-A Simplified Method for Estimating the 
Performance of Supersonic Ram-Jets 


A. B. P. BEETON 


(National Gas Turbine Establishment) 


SuMMaRY: Curves are presented which, together with a knowledge of the intake 
performance, enable ram-jet thrust coefficients and specific consumption figures 
to be obtained quickly and accurately in the range M = 1-0 to 3-0. 


1. Introduction 


The component parts of a typical ram-jet are indicated diagrammatically in 
Fig. 1. For the purposes of performance estimation, it is convenient to consider 
the cycle in two stages, the division being made at Station 1. 


The first phase therefore consists of the air intake, in which the flight kinetic 
energy is converted into pressure. The total temperature 7,, at the diffuser exit 
can be calculated from flight conditions only and does not depend on intake per- 
formance. The corresponding total pressure P,,, however, depends also on the 
aerodynamic design of the intake and on the mass flow demand, and therefore 
involves considerations which come outside the scope of the present paper. It has, 
therefore, to be assumed that pressure recovery characteristics are available else- 
where for the intake in question. 


In the remainder of the ram-jet cycle, the air is first passed through the 
stabilising baffles, thereby suffering a certain loss of pressure. Then it passes 
through the parailel combustion chamber in which a typical hydrocarbon fuel is 
burnt, and in which there is assumed to be no further friction drag on the walls. 
The resulting gases are finally expanded isentropically through a convergent- 
divergent nozzle. 


As the ram-jet is fundamentally so very simple a form of propulsion system, 
it might be expected that the presentation of its thrust characteristics would be an 
equally simple matter. Perhaps the main reason why this is not so is that the high 
gas temperatures used mean that the concept of a constant specific heat (and at 
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Ficure 1. Diagram showing the components of a ram-jet engine and the notation used. 


the very highest temperatures indeed the concept of a constant gas composition 
also), in terms of which the performance of other types of power plant can be so 
successfully expressed, is no longer valid. 


In attempting to determine the theoretical thrust direct from fundamental 
data, it would therefore first be necessary to make a complicated analysis of the 
combustion gas composition at various parts of the cycle, and then to calculate 
therefrom the mean molecular weights, pressures, temperatures, velocities, and so 
on. The majority of these quantities, however, are of no direct importance from 
the performance point of view. A very great saving in time is therefore to be made 
by a graphical presentation which eliminates all or most of these superfluous 
variables. The aim is to devise a system using the minimum number of parameters 
which can both be read off curves against the basic ram-jet design data (area 
ratios, baffle loss, and so on), and which will also enable the required thrust 
coefficient to be calculated from a simple formula. 


In addition, the presentation should be such that there is no significant loss 
in accuracy due to any interpolation necessary, or by having an unduly high ratio 
of maximum ordinate to minimum ordinate in any one diagram. A further desirable 
object is thought to be a treatment that is as far as possible non-dimensional in 
application, so that potential users do not have to cast aside their favourite systems 
of units. 


NOTATION 
A_ flow area 


C, thrust coefficient 


mv 


3 +AP 


stream thrust, = 


standard gravitational acceleration 


air mass flow per second 
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Mach number 

static pressure 

total pressure 

gas constant for air 
total temperature (°K) 
velocity 

fuel/air ratio 

ratio of specific heats for air, = 1-400 
combustion efficiency 
gross thrust 

net thrust 

baffle loss coefficient 
density 


Suffixes 
These correspond to the stations shown in Fig. 1. 


2. Definition of Stream Thrust F 


For ram-jet calculations it is convenient to introduce a flow function F known 
as the stream thrust and defined as 


pa « « « 


This function has the property of remaining constant during frictionless flow in 
a duct of constant area, however much heat is added. With the suffix notation 
corresponding to the stations shown in Fig. 1, therefore, 


F,=F,. 


A further useful property is that the value of F goes through a minimum during 
isentropic expansion through a choked nozzle, and this minimum value can be 
predicted by theory as a function of the air inlet temperature and the fuel/air ratio 
in the combustion chamber. Since, moreover, it is only necessary to measure the 
test bed thrust developed under the prescribed conditions in order to deduce the 
actual value of F at the nozzle throat, a convenient yardstick of combustion 
chamber performance is thereby available. 
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3. Formula for the Thrust Coefficient 


With reference to an atmospheric datum pressure P,, the gross thrust at the 
exit section is given by 


+A, (P,—P,). 


m(1+2)v, 


Since the value of the stream thrust F at this point is by definition 


+A,P,, 


The net thrust 6, is the difference between #, and the momentum of the intake 
air under free-stream conditions, i.e. 


The thrust coefficient C; is usually evaluated on the maximum area A,, so 
that 


Substituting the expression for 6, from equation (3), and putting y= 1-400, this will 
be found to reduce to 


me P.) Gea) 


4. Formula for the Nozzle Exit Pressure 


The expression for Cy given in equation (5) was derived on the assumption 
that the nozzle exit pressure P, was either greater than the atmospheric ambient 
pressure P,, or at least not sufficiently below it to cause the supersonic flow in the 
diverging part of the nozzle to break down. It will also have been observed that 
P, is not evaluated in the process of calculating Cy, and so a check may have to be 
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made periodically to ensure that the jet is not being over-expanded. It is in fact 
quite easily found from the relation 


SIE) 


Since AxP,,/Fx=0-80+41 per cent for all air inlet temperatures and fuel/air ratios, 
and there is certainly no need to work out P, to any greater accuracy than this, we 


Formula for the Specific Fuel Consumption 


The specific fuel consumption (per hour) is given by 


3,600am _ —_7,200am 


Derivation of the Curves 


It is not proposed to explain here in great detail the methods by which the 
combustion data of Figs. 2, 3, 4 and 5 were derived. The following brief description 
indicates the general procedure. 


The fuel was assumed to consist of 15 per cent hydrogen and 85 per cent 
carbon, with a lower calorific value of 10,350 C.H.U./Ib.; this can be taken as typical 
of a liquid fuel of the kerosine type. Combustion was supposed to occur in an 
atmosphere made up of 76-8 per cent nitrogen and 23-2 per cent oxygen by weight. 


The combustion gases were assumed to consist of an equilibrium mixture of 
N., O,, CO,, CO, H,O and H,, for which total heat and entropy data were 
available as functions of temperature every 100 degrees up to 2500°K. Equilibrium 
constants consistent with these values were also worked out over the same range of 
temperatures. At various gas temperatures and pressures for a given fuel/air 
ratio, the first thing to do was to work out the gas composition, and hence the mean 
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Ficure 2. Calculated air specific stream thrusts. 
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FicureE 3. Stream thrust ratios for convergent nozzle. 
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Ficure 4. Stream thrust ratios for divergent nozzle. 
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Ficure 5. Pressure ratios for diverging nozzle. 
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molecular weight, total heat and entropy. This enabled a total heat/entropy 
diagram to be plotted in the region of one atmosphere total pressure and at a 
level of total heat corresponding to air inlet temperatures between 300° and 600°K. 
The next step was to use this diagram to find the variations in gas velocity, 
temperature, pressure and mean molecular weight during isentropic expansion 
from the various input levels of total heat. Finally these values were used to 
deduce the required parameters F/m and A/Ax down the expansion line. A series 
of such expansions was considered at total pressures of 4 atm., 2 atm., 1 atm. 
} atm., and } atm., at air inlet temperatures of 300°, 400°, 500° and 600°K, and 
at eight fuel/air ratios between 0-02 and 0-07. 


Any measurable effect of absolute pressure on the flow parameters is limited 
to those fuel/air ratios close to chemically-correct in which dissociation is both 
considerable and also varies as expansion proceeds. It was found that a maximum 
error of about 4 per cent resulted from using the 1 atm. curve at either 4 atm. or 
4 atm., which represent the approximate practical limits of operation. For 
fuel/air ratios less than 0-05, the error was negligible. 


Figure 2 shows the minimum (i.e. throat) value of F/m down the expansion 
line as a function of 7,, and the theoretical fuel/air ratio «7 at a total pressure of 
one atmosphere. 


Figure 3 gives the required data for expansion in the converging portion of 
the nozzle as a single curve of F,/Fx against A,/Ax. It was found that this would 
cover all the fuel/air ratios and air inlet temperatures investigated +1 per cent, 
which was considered good enough for the present purpose. The maximum area 
ratio considered was 2°4 to 1. 


Similarly Figures 4 and 5 show, for the diverging part of the nozzle, the 
corresponding ratio F,/Fx and the total to static pressure ratio P,./P,. The same 
range of area ratios was covered as for the converging part of the nozzle. 


It will always be possible to sacrifice accuracy in order to gain in simplicity 
and, in the present instance, the unique curves given in Figs. 2, 3, 4 and 5 could 
have been elaborated upon to give corrections for different fuel/air ratios, air inlet 
temperatures and absolute pressures. As already stated, exact values were worked 
out over the range of «;=0-02 to 0:07, T,,=300°K to 600°K and P,,.=4 atm. to 
} atm. abs. The maximum differences between these theoretically exact values 
and the curves in Fig. 2 were found to be +4 per cent (due to the effect of pressure 
only in this case), +1 per cent in Figs 3 and 4, and +14 per cent in Fig. 5, which 
is only used as a relatively rough check. 


It would be unrealistic to suppose that small effects such as these have any 
significant importance in performance estimations, and they have accordingly been 
suppressed here in view of the very considerable gain in simplicity of presentation 
that is thereby obtained. 


Figure 6 merely expresses the baffle loss in terms of quantities which can be 
used directly in the thrust coefficient formula: it does not depend in any way on 
the thermodynamics of combustion. 
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FiGureE 6. Baffle loss curves. 
Defining A as the pressure loss per unit inlet velocity head, we have 
r 2 


Since the mass flow, total temperature and flow area are all unaltered, it is possible 
to write A in terms of the Mach numbers M, and M, in front of and behind the 


baffles respectively: — 
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The stream thrust after the baffles is given by 


mv 
+A,P., dyna 
man} 
which will be found to reduce to requi 
reco\ 
knov 
praci 

Finally, in front of the baffles a flow function can be evaluated in terms of a 

M. as IStic. 
agai 
AP, 1 opel 
Combining this with the previous expression gives 
F,  M, (1+yM,?) 
Now it is known that the value of F after the baffles is equal to the value found crit 
for the nozzle at the other end of the combustion chamber. Hence calculating pre 
coe: 
mi (Ty) mV 
specifies M,, from equation (10). So if the baffle loss A is known, M, is uniquely ¢ 
determined by equation (9). Finally with both M, and M, known, the value of 
A,P,./F,;=A,P,:/F, is given by equation (12). 

Since M, and M, have no significance from the performance point of view, they It 
can be eliminated and the results merely expressed as values of A,P,./F, for A, 
different values of A and F,/{m4/(T,,)}, as in Fig. 6. 

The quantity (1/M,*)(P,./P,) given in Fig. 7 is a simple function of Mach 
number. It is a more suitable function to plot than P,./P, in that it covers a 
relatively small range of values (between 1:6 and 4:1) for M, between 1 and 3, 
and can therefore be read off the curve reasonably accurately. Otherwise it can 
be evaluated to any accuracy desired from the relation 

y/G-1) wl 


p= (1+40-DM") 


Similarly v,/{g/(T,:)} is plotted as another function of M, in Fig. 8:— 
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1. Notes on Using the Curves 


In general, the pressure recovery factor P,;,/P,, is determined by the aero- 
dynamics of the intake, as affected by the degree of matching with the engine. For 
many purposes the intake can be considered to be ideally matched to engine 
requirements at all the points it is necessary to investigate, and the design-point 
recovery can then be expressed straight away as some function of the flight Mach 
number. The likely variation of this maximum recovery with M, is fairly well 
known, but typical values would be given with more authority in a book on 
practical aerodynamics than in a purely theoretical treatment such as this. 


When off-design conditions are required with an intake of known character- 
istics, it is best to re-plot them, if necessary, to show the flow function A,P,./m 
against P,,/P,.. The value of this flow function can be found from the engine 
operating conditions as the product 


which thus enables P,,/P,, to be picked off the intake characteristic. 


Another usual case is one in which the intake is known to be running super- 
critically and the critical condition at the given flight Mach number has been 
previously worked out, so that (A,P,:/m).i: can be determined. The thrust 
coefficient then becomes 


Pe F,\ (Fx\ (T.) (2 1 A, ax 


It is therefore unnecessary to work out either P,./P,., F,/Fx, F;/{mV(T\:)} or 
A,P,,./F, in order to derive Cy in this case. 


Similarly the specific fuel consumption formula (7) reduces to 


AP (14) 
while the expression for the nozzle exit pressure is now 
Po\ (Fx 


(Abe) (Fe) 
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Figure 2 shows the ratio Fx/m, also known as the air specific stream thrust, 
in terms of the air inlet temperature T,, and the effective or theoretical fuel/air 
ratio an. It will be seen how this parameter approaches a maximum at approxi- 
mately chemically-correct mixture. In practice, however, the greatest value 
measured on the test bed is somewhat less than the theoretical limit at the given 
air inlet temperature, and the peak performance of the engine can conveniently be 
expressed as the maximum air specific stream thrust actually developed. Calcula- 
tion of Cymax then proceeds directly from this value of Fy/m, without reference 
to Fig. 2. 


In general, it is clear that the combustion chamber operating conditions can be 
specified by an one of the following : — 


(i) Fx/m and a, 
(ii) Fx/m and », 
(iii) and 2, 
(iv) ay and », 
(v) @ and ». 


For calculating the thrust coefficient, it is only necessary either to know Fx/m or 
to be able to obtain it from a given value of ay from the curves in Fig. 2. For 
calculating the specific fuel consumption, the value of the true fuel/air ratio « also 
has to be deduced from the further information supplied in any of the five ways 
listed. 


It should be noticed, however, that those quantities obtained directly from 
test-bed results, namely Fy/m and z, are also those which lead to the most direct 
determination of thrust coefficient and specific consumption. It would neverthe- 
less be more reasonable to specify instead some typical value of 7 when considering 
a “hypothetical ” engine of unknown characteristics. 


Of the remaining terms in the expressions derived for Cy, P,/P, and specific 
consumption, it will be seen that F,/Fx, F,/Fx, Ps:/P,, (1/M,)(Po/P,), and 
vo/{gV(T,:)} are obtainable directly from Figs. 3, 4, 5, 7 and 8 respectively. 
F,/{m+/(T,,)} then has to be calculated from the three quantities T,., F,/Fx and 
Fx/m. As well as appearing in the formula for C;, this parameter is used to evaluate 
A,P,./F;, the only quantity so far undetermined, from the curves in Fig. 6. 


It will be observed that the only units specified in this paper are the second 
and the degree Centigrade absolute. All the figures and formulae are therefore 
equally well suited for working in metric as in ft.-lb.-sec. units. 


8. Example 


Any doubts as to the exact method of using Figs. 2 to 8 should be resolved 
by following through the example in this section. 
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Data 


Flight Mach number M,=2:30 
Inlet total temperature T .=450°K 
Intake pressure recovery ratio Py. / Po =0-75 
Baffle loss coefficient A=6 
Effective fuel/air ratio an=0-045 
Combustion efficiency n=0-80 
Nozzle contraction ratio A,/Ax=1-50 


Nozzle divergence ratio A,/Ax= 1-30. 


Hence : — 


from Fig. 2, Fx/m=150-0 sec. 
from Fig. 7, (1/M,”) (P../P,)=2°36 
from Fig. 8, vol {gv (T =3-275 sec. (°K)-?/? 
from Fig. F,/Fx=1-322 
from Fig. F,/Fx=1-078 
and from Fig. P;./P,= 4-20. 


Therefore 


Hence from Fig. 6, 


Therefore Cy= > 


0:75 x 2:36 { 1:078 oar} 1:30 


Loris 934 J ~ 
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“Co 


_ 0°800 x 2-30 x 2-30 x 0°75 x 2:36 x 1:50 
4:20 x 10775 x 1-322 


= 1°88, 


indicating that there is still plenty of pressure ratio in hand at the nozzle exit. 


Finally, the specific fuel consumption is 


5.143 


5,143 x 0-045 x 2:36 x 0-75 
~ 0-80 x 0°857 x 150-0 x 1-322 x 1-0775 


= 2:80 Ib./hr./Ib. (or kg./hr./kg.) 
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Variably Reinforced Circular Holes in 
Stressed Plates 


RAYMOND HICKS, Ph.D., M.Sc. 


(G.E.C. Simon-Carves Atomic Energy Group) 


Summary: Expressions are obtained for the radial, tangential, and shear stresses 
in an infinite plate containing a variably reinforced circular hole. The variation 
in the cross-sectional area of the reinforcement is chosen so that the maximum 
tensile resistance of the reinforcement occurs at points around the hole where 
the stresses in the plate tend to be high. It is shown that when the reinforce- 
ment has a given weight, the effect of varying its cross section is to reduce the 
stress concentration in the plate. 


The following particular applied stress systems are considered : — 
(i) Uniform tension or compression in one direction. 


(ii) Uniform shear. 


(iii) Unequal principal stresses. 


For the particular cases of uniform end loading and unequal applied stresses 

in the ratio 2: 1, data have been obtained which enable the designer to 

determine the stress distribution in plates with reinforcements having a 
reasonably wide range of practical dimensions. 


1. Introduction 


The presence of a circular hole in a flat plate loaded in its own plane results in 
a stress concentration. If the hole is small compared with the lateral dimensions of 
the plate its effect on the stress distribution will be local, whereas a large hole may 
appreciably alter the stress distribution in the whole plate. For example, Kirsch” 
found that the effect of a circular hole in an infinite plate uniformly stressed in one 
direction is to cause a stress concentration of three, while Timoshenko™ showed 
that Kirsch’s solution can be applied to a finite plate provided that its least lateral 
dimension is several times the diameter of the hole. Thus, when the width of the 
plate is four times the diameter of the hole, the error in the solution when calculat- 
ing the maximum stress is six per cent. 


A bibliography of some important papers dealing with unreinforced circular 
holes in stressed plates is given in Ref. 3. 
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It has been shown by Gurney“, Timoshenko’, Beskin“’ and many others (sce, 
for example, Ref. 3) that the effect of suitably reinforcing the hole is to decrease 
the maximum stress in the plate. Thus, when the maximum stress is nowhere to 
exceed a given value, the presence of the reinforcement results in a decrease in the 
necessary weight of the structure. The design of the reinforcement is therefore an 
important problem when the weight of the structure must be as low as possible, for 
example, with holes in the skin of an aircraft. 


It is worth noting that Gurney used a relatively thin annular type of reinforce- 
ment and determined the stress distributions in the plate and reinforcement due to 
various load systems applied at infinity. Timoshenko, however, used a compact 
ring as a reinforcement, the radial width of which was small compared with the 
diameter of the hole. That is, a reinforcement with a negligible bending stiffness in 
the plane of the plate. 


Reinforced holes which do not alter the stress distribution in a plate are said 
to be neutral. The problem of a neutral hole in a plate loaded in its own plane has 
been considered by Mansfield“ who found that such holes are possible if there is 
a choice for the shape of the hole and the tensile stiffness of the reinforcement. 


For many practical problems the necessary cross-sectional area of the rein- 
forcement which makes the hole neutral is sufficiently large to be prohibitive, while 
the necessary shape of the hole in some cases is defined by a curve which does 
not close. 


Mansfield also considered the problem of a variably reinforced circular hole 
in an end-loaded plate. Using an inverse method of solution and taking stress 
functions which gave two independent stress pertubations around the hole, he found 
expressions for the cross-sectional areas of reinforcements corresponding to the 
chosen number of stress pertubations. He found that the necessary cross-sectional 
areas of the reinforcements were defined by 


Ag = cos nd 


where n is 0, 2 or 4, depending on the applied stress distribution. 


The problem of a reinforced elliptical hole in an end-loaded plate has been 
considered by Wells“ for several important applied stress distributions, while the 
author’ has considered the design of similar reinforced holes for stress systems 
of the type found in cylindrical and ellipsoidal pressure vessels. 


This investigation deals with a reinforced hole in an infinite plate under the 
action of the following applied stress systems : — 


(i) Uniform tension in one direction. 
(ii) Uniform shear. 


(iii) Unequal principal stresses. 
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The reinforcement has a variable cross-sectional area, so designed that its tensile 
resistance increases at points around the hole where the stresses in the plate tend 
to be high. The advantage of a reinforcement of this type is that no material is 
wasted in reinforcing parts of the hole where the accompanying stresses are 
relatively small. 


Expressions are obtained for the radial. tangential, and shear stresses in the 
plate. These expressions are written in terms of a set of constants which can be 
evaluated readily for any particular numerical problem. As would be expected, 
the concentration of stress is local. so that the solution can be applied to a finite 
plate provided that its least lateral dimension is several times the diameter of 
the hole. 


For the problem of a plate uniformly pulled in one direction, the constants 
appearing in the expressions for the stress components have been evaluated for a 
wide range of reinforcements with practical dimensions, it being assumed that 
Poisson’s ratio is 0-3. Furthermore, curves have been obtained for the maximum 
tangential stresses around holes in plates having reinforcements of a constant 
weight, but with different variations in their cross-sectional areas. It has been 
found that, for a given maximum tangential stress, the weight of the reinforcement 
can be reduced by as much as 50 per cent by suitably varying its cross-sectional 
area. 


Similarly, for the important case of unequal applied principal stresses in the 
ratio 2:1, values of the constants appearing in the expressions for the stress com- 
ponents have been calculated for the same range of practical reinforcements. 


NOTATION 


stress components 


f..fy principal stresses 
applied stress 
applied shear stress 
tensile force in reinforcement 
polar co-ordinates 
Cartesian co-ordinates 
radius of hole 
thickness of plate 
area of variable reinforcement 
area of constant reinforcement 
stress function 
Young’s modulus 
M,N _ parameters defined in equations (19) and (36) 
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Poisson’s ratio 
strain 
a’, B’,2,8 constants 


R. parameters introduced in Section 7. 


Suffix 


b denotes boundary conditions. 


2. Uniform Tension 


For an infinite plate uniformly stressed in the direction of the x-axis, the 
boundary conditions at infinity are 


(rr). = + cos 28) | 


sin 20 


where the @ co-ordinate is measured from the direction of the x-axis. 


Referring to Fig. 1, around the hole the boundary conditions are obtained by 
considering an elemental length of the reinforcement subtended by an angle dt. 
There is assumed to be no radial shear on a radial cross section of the reinforce- 
ment, because its berding stiffness in the plane of the plate is small. This follows 
from the work of Reissner and Morduchow”, who concluded that reinforcements 
which have a relatively small bending stiffness in the plane of the plate experience 
primarily tensile and compressive stresses as opposed to bending stresses. 


a8 
FiGure |, Forces acting on an elemental length of reinforcement. 
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For force equilibrium in a radial and tangential direction respectively, it can 
be shown that 


F=bh (rr), 


— hb (r6), | 


Also, equating the strain in the reinforcement to the tangential strain in the adjacent 
part of the plate, it is found that 
F 
= (06 — rr), (3) 


R 


it being assumed that the materials of the reinforcement and the plate are the same. 


Combining (2) and (3), the following boundary condition is obtained 
va), 22 (4) 
b bh b . 


From Ref. 1, the stresses around an unreinforced hole in a plate uniformly 
stressed in one direction are 


(rr), =0 | 


(r6),=0 
(46), =f (1 —2 cos 26) | 


Thus the presence of the hole increases the maximum stress in the plate by a factor 
of three. The effect of suitably reinforcing the hole is to decrease the stress 
concentration in the plate, and in this analysis the reinforcement is designed so 
that its cross-sectional area varies in a similar manner to the tangential stress given 
by (5), that is 


where 2’ and f’ are constants so that z’+ 6’ is the area of the reinforcement at the 
point 6=0, and 2’—/ is the area of the reinforcement at the point @=7/2. It is 
clear that z’ is the cross-sectional area of a reinforcement having the same weight 
as the variable reinforcement defined by expression (6). Comparing (5) and (6), it 
can be seen that the advantage of the variable reinforcement is that its effectiveness 
in decreasing the tangential stress around the hole increases as the tangential stress 
tends to increase. 


Using the expression for A, given by (6), the boundary condition specified by 
(4) can be written as 


where z has been written for 2’/bh and 8 for £’/ bh. 
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If # is a function which satisfies the equation 


the stresses in the plate are given by 


~ 


0 
06 


The solution of (8) which can be made to satisfy the boundary conditions is 


n=2 


where n is even and A,, A,’, B,’, etc., are arbitrary constants. It is worth 
noting that if A,, B,, C, and D, are zero and sin n@ is substituted for cos n@, the 
resulting stress function can be used to solve the problem of a variably reinforced 
hole in a plate having a pure shear distribution applied at infinity. 


Substituting (10) in (9), the stress components become 


+B,(1+2logr)+2C,-— [ n (n—1) A,’r"~2+(n +1) (n— 2) B,’r" + 
A, 


+tn(n+l) 


+(n+-2)(n—1) | cos né 
A x 


+(n+2)(n+1) B,’r"+n(n+ 1) +(n—2)(n-1) =| cos né 


ro = n| (n- 1) A,/r"-24+(n+ 1) 1) ~(n-1) sin 


Using these expressions, for the stresses to remain finite as r approaches infinity it 
is found that 


B,=0 
A,/=0 when n+ 2 (12) 


B,’=0 when n>2 
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REINFORCED HOLES IN STRESSED PLATES 


Also, from (1), ‘ (13) 


Eliminating F from the equilibrium equations given by (2), the following boundary 
condition is obtained 


This condition is satisfied when 


fe, An 


where 6.=1 when n=2 


5,=0 when 


Using the values of the constants given by (12), (13) and (15), the stress components 
become 


n=2 


{ A, 1) (n- fb. cos no 


J 
2 


It is worth noting that these expressions are independent of the cross- 
sectional area of the reinforcement. They are, therefore, valid for any type of 
reinforcement provided its bending stiffness in the plane of the plate is negligible. 
Furthermore, by applying the principle of superposition, these expressions can be 
used to find corresponding expressions for the stress distribution in a plate 
subjected to a pure shear or an unequal principal stress distribution at infinity. 


The arbitrary constants appearing in the expressions for the stress components 
must be chosen to satisfy the boundary condition given by (7). Thus, substituting 
the appropriate expressions for (rr), and (68), in (7) and equating to zero the 
coefficients of each cos n@ term in the resulting form of the boundary condition, 
the following set of equations is obtained for the arbitrary constants A, and A,. 
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Aue] <0} 
2. 4, 


where N= ty) 24 1]- 


Bf A, f 
M= (l—v)— B(1+yv)— 2-1) 


when n= m, when n=m, 
when A_.=0. 


Yn=1 when n~™m, 


From (4) and (16), it follows that, for any given set of values for the constants A,. 
the solution of the problem is exact provided that the reinforcement is chosen so 
that the variation in its cross-sectional area is given by the expression 


Ar v) sin?é — (1 +2 3 A+ cos ne } 


n= 2 


=bh { +fsin*6—2 cos no} (20) 


Thus, if there is no restriction on the variation of the cross-sectional area of the 
reinforcement, and if there is only one particular stress system applied to the plate. 
the expression (20) can be used to design variable reinforcements corresponding to 
any reasonable number of stress perturbations. However, when dealing with 
particular problems the constants A, must be chosen to satisfy certain design con- 
ditions. For example, the chosen values for the constants must correspond to 
reinforcements having practical dimensions and with stress concentrations which are 
relatively low for a reinforcement of a given weight. As previously stated, this 
method of solution has been used by Mansfield. 


In this analysis it is proposed to use an inverse method of solution. Thus, 
when dealing with practical problems it is assumed initially that the stress distribu- 
tion in the plate is given with sufficient accuracy by taking a finite number of terms 
in the expression for the stress function. The constants A, are then determined 
from (18) and the cross-sectional area of the reinforcement necessary to make the 
solution exact is determined from (20). Finally, it is verified that the cross-sectional 
area of the reinforcement calculated from (20) does not differ appreciably from the 
given cross-sectional area of the reinforcement defined by (6). 
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REINFORCED HOLES IN STRESSED PLATES 


As an example, consider the following extreme case where @=—a and 
2=0-1. Assuming A,=0 when n> 8, the solution of the first five of equations 
(18) gives 


o= —0:352182 b*f 
»= —0°150104 b*f 
, = —0-020629 b°f 
—0-003843 b*f 
s= —0:000858 b''f 


A,= 
A, 
A 
A,= 
A,= 


Substituting from (21) into (20), the stress distribution associated with these 
constants is exact when the variation in the cross-sectional area of the reinforcement 
is given by (20). It can be verified that this distribution in cross-sectional area of 
the reinforcement is virtually the same as 


bh (1 —cos 28), 


which is the given expression for the variation in the cross-sectional area of the 
reinforcement. Thus, for this particular example. the inverse method of obtaining 
the arbitrary constants A,, is justified. 


Substituting from (21) into (16), it is found that, when r is several times greater 
than the hole diameter, the stress components reduce to 


= (1+ cos 24) 


f 


(1—cos 24) 


9 
5 Sin 26. 


That is. the effect of the hole is local, and for this reason the solution can be applied 
to a finite plate provided that its width is several times greater than the diameter 
of the hole. 


3. Design Data for Uniform Tension 


For design purposes, values of a, (a, f b"*?= —A,) are given in Table I for a 
wide range of reinforcements with practical dimensions. Using these values for the 
constants, it can be verified that the cross-sectional areas of the variable reinforce- 
ments defined by (20) are virtually the same as the cross-sectional areas of the 
reinforcements defined by (6). That is, the inverse method of determining the stress 
distribution is justified for all reinforcements given in Table I. 
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TABLE I 


VALUES OF dn (a, b" +? f= — An) FOR VARIOUS VALUES OF 7 AND 8 WHEN v=0°3 


—0-7 


():35218 0°36321 036883 0°38035 
0:15010 0° 15368 0°15574 016036 
0:02063 0:01639 0:01435 0:01032 
0:00384 0:00237 0:00179 0-00090 
000086 000042 0:00028 0:00010 


0:19992 0°21485 0:22280 0:23776 0:25470 
0:06526 0:06864 0:07112 0:07729 0°08492 
0:02179 0:01668 0:01451 0:01047 0:00650 
0:00646 0:00370 0:00275 000136 0:00049 
0°00196 0:00090 0:00058 0:00020 0:00004 


0°11796 0:13104 0:13745 0°15122 0°16715 
0:02598 0:02938 0:03200 0:03846 0:04638 
0:01764 0:01331 0:01160 0:00849 0:00537 
0:00596 000332 0:00246 0:00122 0:00045 
0:00195 0:00087 0:00056 0:00020 0:00004 


0:04558 0:05359 0:05782 0:06793 0:08097 
~0°00492 0-00062 0:00220 0:00865 0:01619 
0-00226 0:00943 0:00833 0-00627 0:00408 
0:00450 0:00253 0:00189 0:00096 0:00036 
0:00154 000069 000045 0:00016 0:00004 


It is worth noting that when z= £=0 the hole is unreinforced and the solution 
of (18) gives 


fb* 
4 


A,=0 when n>2 


A,=— 


Also, for a reinforcement with a constant cross-sectional area, 8=0 and the 
expressions for the constants are 


A,=0  when..2n> 2. 
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Thus, around the hole 
af 7 22f cos 26 
{l+2(1+yv)} [l+2(3+y)] 


oy, 2(1+2»)f cos 20 
[1+2(1+)] [l1+2(3+9)] 

— 4af sin 26 

[l+2(3+%)] 


(rr), = 


= 


The expressions show that, when there is no reinforcement, the maximum stress 
acts in a tangential direction on the boundary of the hole and equals 3f. Further- 
more, this maximum stress decreases as the area of reinforcement increases from 
zero. 


It is outside the scope of this paper to determine the maximum principal 
stresses or octahedral stress concentration factors for a wide practical range of 


MAXIMUM TANGENTIAL STRESS 
=h 


FiGure 2. Maximum tangential stress for constant weight reinforcements. 
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reinforcements. Furthermore, for any particular example these can readily be 
obtained by using the information given in Table I. However, to illustrate the 
effect of varying the cross-sectional area of the reinforcement, curves showing the 
variation in the maximum tangential stress for constant weight reinforcements are 
given in Fig. 2. From these curves it can be concluded that for a given maximum 
tangential stress, the necessary weight of the reinforcement can be reduced by as 
much as 50 per cent by suitably varying the cross-sectional area of the 
reinforcement. 


4. Pure Shear 


Using the principle of superposition, and rotating the co-ordinate axes through 
45°, the stress components g’ven by (16) can be used to obtain the following 
expressions for the stress components in an infinite plate under the action of a 
constant shear stress distribution. 


Thus, if g is the applied shear stress 


{ [ n(n+1) —(n- | } sin né 


A 
[n 2)| } sin n@ 


x A, q hb? } 
N 
[ (n+ -& +4 (1 + cos né. 


The stresses around an unreinforced hole in a similarly loaded plate are 


(rr), = =0 


(27) 
(60), = — 4q sin 26 


Thus the stress concentrations occur on radial lines which bisect the co-ordinate 
axis. It therefore follows that a suitable reinforcement would have a variable 
cross-sectional area given by 


From (4) and (28), the compatibility condition becomes 


(60 — vr), (2+ B cos 46)—(rr),=0, . 
where z and / are given by z’/bh and ’/bh respectively. 


Substituting in (29) the appropriate expressions for the stress components, and 
equating to zero the coefficients of the resulting sinm# terms, we obtain the 
boundary condition as given by equations (30). 
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— 


n-4 n+4 
3 (2n 941) + BAH of 0, 2.4 


where b6n=1 when m= an, 6,=0 when ma, 


A_,=A_.=A,=0. 


As previously shown, by using an inverse method of solution the set of equations 
(30) can be used to determine the constants A, for any particular problem. Also, 
since the stress distribution is symmetrical with respect to lines inclined at 45° to 
the directions of the applied shear, the constants A,, A,, A,... . . must be zero. 


From (4) and (26), it follows that, for the values of A, obtained by using an 
inverse method of solution, the solution of the problem is exact provided that the 
variation in the cross-sectional area of the reinforcement is given by the expression 


bh { qsin 26+2 nt} 
{ sin 26 1+} 


Finally, for a reinforcement with a constant cross-sectional area, 8=0 and the 
constants are given by 


A,=0 when n+2. 


5. Unequal Principal Stresses 


Again using the principle of superposition, it is found from (16) that, for a 
plate under the action of eel stresses f, and f, at infinity, the stress components 
are given by 


fatfy {4 [mins + 


~ ; } cos né 
{45 [mins —(n-1)(n- 2] - 
_ 


cos né 
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The stresses around an unreinforced hole in a similarly loaded plate are 


(rr), = (78), =0 


(34) 
(60), =(f,+f,) —2 (f.—f,) cos 24 


A suitable reinforcement would therefore have a cross-sectional area defined by (6), 
in which case the compatibility condition given by (4) can be reduced to the follow- 
ing set of equations : — 


{ [2 Qn-1+9)+1]7, +B[ Qn-54) 43 7, - 


8, 


fy) 
2 


M= 


and y,=0 when m=n, when ma, 


A..=0. 


From (4) and (33), it is found that the solution of the problem is exact when the 
cross-sectional area of the reinforcement is given by 


—v) cos 26+2 = [2n— 1+ ¥] cos nt | 


=bh | x cos nt (37) 
n=2 


When the reinforcement has a constant cross-sectional area, the solution of (35) 
gives 


__ 


4 


and A,=0 when n > 2. 
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TABLE II 


VALUES OF @n (@nb"+?fy=An) FOR VARIOUS VALUES OF x AND 8 WHEN v=0°3 AND fr=2fy. 


—0°5 

—1:19548 —1:19899 —1-:20155 — 120826 
—0-08064 -—0:09900  —0-°10820 0°12674 
—0:01357. -—0:01206 -—0-01108 — 000870 
—0°00253 -—0:00174 —0-00138 — 0:00076 
—0:00056 -—0:00031 —0-00022 — 0:00084 
—0°86146 —0°83843  —0°83173 —0°82621 
0:06559 002830 001133 —0-02083 
0:00231 —0-:00339 —0-00487 —0:02593 
0:00069 -—0:00075 -—0-00092 —0-00077 
0:00021 —0:00018 —0-00020 —0-00012 


—0°66992 —0°61289 —0-59549 — 0°57566 —0°57162 
0: 13204 0:08050 0:05957 0:02253 —0-01129 
0:01764 0:00444 0:00089 —0:00283 —0-00348 
0:00596 000111 0:00019 —0:00041 —0-00029 
0-00195 0:00029 0:00004 —0:00007 — 0:00003 


—0°48956 —0°38573  —0°35634 ~0°30960 
0°18133 0°11311 0:08925 0°01716 
0:03269 0:01081 0-00547 —0-00211 
0:01201 0:00290 0:00124 —0-00019 
0:00410 0:00079 0:00030 —0-00002 


Design Data for f,=2f, 


For the important practical case when f,=2f,, values of a, (a,f,b"*?=A,) are 
given in Table II for a wide range of reinforcements with practical dimensions. As 
in the problem of a uniform tension in one direction, these values for A, have been 
found by using an inverse method of solution. 


1. Parameter Solution for Unequal Principal Stresses 


When the plate is under the action of a principal stress system of the type 
found in cylindrical and ellipsoidal pressure vessels, reinforcements having variable 
cross-sectional areas, which result in reductions in the stress concentrations, can be 
designed by using the parameter method described in a previous paper” dealing 
with reinforced elliptical holes. 


Referring to Fig. 3, the applied stress system shown in State 1 can be resolved 
with the aid of the parameter into the component stress systems shown in State 2 
and State 3. Assuming that for State 2 the plate has a circular hole with a rein- 
forcement having a constant cross-sectional area Acx, the stress components are 
given by equations (39). 
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STATE | 


if 


| STATE 3 


STATE 2 


FicurE 3. Resolution of the applied stress system into two component applied stress systems. 


2 


The suffix 2 indicates that the stress components are associated with the applied 
loading shown in State 2. 


Around the hole, the tangential strain in the plate is given by 


_ (n-o+)Df — 2(n-0-I1)f 
[14+6+y)R] 
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Assuming that for State 3 the plate has an unreinforced hole, the stress components 
become 


rr,= (1- =) cos 26 


Thus, around the hole, 


¢,E=wf (1 —2 cos 26). ‘ ‘ (42) 


Finally, assuming that for State | the stress distribution in the plate is the sum of 
the stress distributions for States 2 and 3, the problem reduces to finding the cross- 
sectional area of the reinforcement necessary to satisfy the boundary conditions 
around the hole. 


From (40) and (42). around the hole the tangential strain in the plate becomes 


[n+1+(1 oR] 
{[1+(1+y) R] 


[n—1+(3+%) oR] 


(1+(3+%) R] 


f-2 f cos 284. (43) 


From (2), the normal force acting on a radial cross section of the reinforcement is 
given by 


F=bh(rr,), 


= { [n-o+I1]R [n—w—-1]R 


From (43), (44), and (45), the variable cross-sectional area of the reinforcement is 
defined by 


Ax{[n+1+(1+ 4) oR] [1 +(3+%) +») R] cos 26} 

= {(n—o + 1) [14+(3 RJ-2 (n-o- 1) [1+(1+%) R] cos 26}. 
Thus the stress distribution in the plate and the variation in the cross-sectional area 
of the reinforcement depend on two arbitrary parameters » and R (note Acg=Rbh). 


This implies that, when dealing with a particular applied stress system, practical 
reinforcements can be obtained by choosing the most suitable values for w and R. 
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TABLE Ill 


VALUES OF A,,/bh FOR VARIOUS VALUES OF R AND © WHEN v=0°3. 


(degrees) | 0-2 


0°184 
0°189 
0°200 
0°208 
0:214 
0°216 
0:217 


0°169 
0°179 
0-200 
0:217 
0°228 
0°235 


0° 156 
0-171 
0:200 
0°225 
0°243 
0°254 


For example. consider the important case of applied principal stresses in the 
ratio two to one, that is f and 2f. Then, if the larger of these two principal stresses 
acts along the radial line where @ is zero, the distribution of tangential stress around 
an unreinforced hole is 


06,={(3-2c0s26. . . . (42) 


Thus, a suitable reinforcement would have a cross-sectional area which was 
maximum at the point @==/2 and minimum at the point 6=0. 


For a reasonably wide range of practical reinforcements, values of Ax/bh for 
various values of R and w are given in Table III. 
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90 0-445 0-682 0-936 1-193 | 
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On the Accuracy of the 
Representation of a Lifting Line 
by a Finite Set of Horseshoe Vortices 


A. R. COLLAR 


(University of Bristol) 


SuMMaARY: The paper examines the error invoived when the continuous distri- 
bution of lift on an unswept lifting surface is represented by a finite number 
of horseshoe vortices, the strengths of which are determined by a prescribed 
downwash. Attention is directed mainly to the case of elliptic loading, but an 
anti-symmetrical case is also studied to ensure that the conclusions are not 
peculiar to the elliptic case. 


If the number of horseshoe vortices is m, it is concluded that 


(a) when the downwash at the centre of each horseshoe vortex is prescribed 
equal to that for the continuous case, the total load due to the horseshoe 
system over-estimates the true load in the ratio (m+ 1)/m, 


(hb) if the circulation is determined for a wing of finite aspect ratio by its 
geometric incidence, the error will be less than that in case (a); a formula 
involving m and the aspect ratio A is derived. 


The analysis involves an unusual method for the reciprocation of the 
matrix of which the typical element is 


1. Introduction 


Many problems in aeronautics arise in which it is convenient to represent a 
lifting surface by a lifting line, and further to represent the lifting line by a finite 
set of horseshoe vortices. Given the strengths of the vortices, the induced down- 
wash is readily written down; but the inverse problem is less easy. 


The present note has been written because there are some points of interest in 
the method of reciprocation of the matrix involved in the solution of the inverse 
problem, and because the final result sheds some light on the accuracy of the 
representation. 


Received March 1958, 
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REPRESENTATION OF LIFT BY HORSESHOE VORTICES 


NOTATION 


A aspect ratio; see also (18) infra 
B, see (27) 

C, see (38) infra 
D, see (39) infra 

triangular matrix 

I unit matrix 

K circulation 

L iift 
Ly rolling moment 
R,, lift ratio 

U square matrix 

Vs fluid speed 

a lift coefficient slope 

c chord; also matrix 
i,j current suffixes 

k_ fraction of width of horseshoe 
m number of horseshoes 
n=m+1 
r,s current suffixes 

semi-span 

square matrix 

w downwash velocity 

@ incidence 


p fluid density 


2. The Problem 


Suppose we have a set of m horseshoe vortices, of strengths K,, K,,.. .. Km, 
the bound vortices, each of length 2s/m, lying on a straight line normal to a stream 
of density p and speed V (Fig. 1). This is perhaps the simplest representation of 
an unswept lifting wing of span 2s. The velocities induced by this system of 
vortices are readily written down; normally we are concerned mainly with the 
downwash on the lifting line itself. Let w; be the downwash at a point on the i 
bound vortex, distant 2sk/m from the left-hand end of that bound vortex, In 
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3. 
I 
K K terms 
|" 
Acco} 
K, K,|-K, Ky|-K, Kr ~Km 
Physi 
DOWNWASH VELOCITIES w, Wy W;----Wm. 
CALCULATED AT POINTS A BC->--M. 
Wn (= Wms) CORRESPONDS TO POINT N. 
FIGURE 1. 
practical cases the fraction k will almost always be given the value 4; but the greater 
generality is retained here. Then we can write at once 
m k+i-1 k+i-2 k+i-m~ k+i-1l—m 
The 
-K —K 
(k+i-1)(k+i-2) (k+i-—m)(k+i-m-1) 
In matrix notation, equation (2) is 
UK, 
m whet 
where the square matrix U. of order m, has the typical element 
(4) 
Given a set of strengths K, equation (3) enables the downwash set w to be 
calculated. The more difficult problem is. however, the solution of Pre- 


K=-—- U-"w, . (5) 


to obtain the set of strengths K for a prescribed downwash set w. We proceed to 
evaluate the matrix U~' where U is defined by (4). 
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3. Reciprocation of U 
If we return to equation (1), we note that w; is expressed as a sum of m+1 
terms, the associated modal column (vector) being 
Kn-Kn-,, 
Accordingly, let us add to the m equations (1) the (m+ 1) equation 


Kn) 
Physically, this defines the downwash at a point just outside the lifting line 
(N, Fig. 1); normally, w,4, will not be prescribed with w. However, we can write 


the full set of equations (1) and (6) as 


8x5 }_ | 1 1 1 || | 
m WwW, | k’ k-1 w k-m K,-0 | (7) 
a 1 | | 
| 1 


The square matrix appearing in (7) is the matrix u of which the reciprocal has been 
given elsewhere by the author’. Since we are now dealing with a square matrix 
of order m+ 1, it is convenient to write 


In this case our symbols accord with those of Ref. 1. The reciprocal of u is v, 
where the typical (r, s)" element is 


(r—1)! (n—-P)! (s—1)!(n—s)! 
Accordingly, from (7) we obtain 


Pre-multiplication by the “triangular” matrix 


0 
| 
ged 
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8xs Theor 
This, while similar to (5), differs in that we have on the right one unspecified 
element w,, and on the left a single specified (zero) element. ] 
Accordingly, write Hv=c expre 
and partition the matrices so that (10) is written ‘ 
K | = Ci Cn || W ; 
| 
| 
0 Cnj Can Wr 
From the last line in (11) we obtain 
and, since c,, is a scalar, the remainder of (11) may be written, on the use of (12). 
K=cyw—- Ww 
Wit 
Cnn 
Comparison of (13) with (5) shows that 
and 
U-'= | ey — (14) 
Can 
It may be remarked that the last row of H may be written [3,5,. . . ., 8), where 
§ is arbitrary: the analysis is not affected. Sul 
4. Evaluation of the Typical Element of U~' 
Although (14) provides a formal solution, and is also a result that may be used 
directly for calculation, it is possible, with a little manipulation, to write down the 
typical element of U~'. To do this, certain theorems, proved in Ref. 1, are wt 
required; they are also required for the further investigation of Section 6. For 
simplicity, the theorems (all particular cases of a more general theorem) are quoted 
here : — 
Theorem B 
TI 
(-ly-' .(k-—r+1) 
Theorem C ( 
rat ly-! 
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Theorem D 


_n(n—1)k(k-1) 


Now, in (14), ci; is the (i, j) element of Hv, which is clearly, in view of the 
expression for H, 


Cy= 


r=1 


G-I)'(n-p! k—r+j 
(15) 


In particular, with i=n the summation in (15) is unity, by Theorem B. so that 


= 


With j=n, the expression (15) becomes 


(k+n—1)....k (kK—r+n). . 


and finally, with j=n, the expression (16) gives 
Can = (n—1)! (18) 


Substitution from (15), (16), (17), (18) in (14) gives, for the typical element of U-’, 


CinCnj _ (k+j-—1). .(k+j—n) 


(-1y-! 1 1 

Thus, finally 
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5. The Special Case k=1/2 


It will be usual, in practice, to define the downwash at the centre of each bound 
vortex, since this preserves symmetry. In this case, k=4, and the matrix U has 
the typical element 


Uu= 27+ 1) 


(20) 


From (19), and in view of (8), the reciprocal will be defined by 


(2j-1). . .(Q2j—2m—-1) (Qm+1-2r). . .(3-2r) 


(21) 


As an illustration, let m=4. Then 


From (21), the inverse is defined by 


(2j-1). . .Qj-9) § @-2r). . .3-2n) 


From this formula, U~' is readily obtained in the form 


| 35.35 15.35 15.21 35.5 | 
35.15 15.95 15.41 35.9 

359 15.41 15.95 35.15 

355 «15.21 15.35 35.35 


It may be noted, for purposes of computation, that there will always be a factor in 
each column corresponding to the terms preceding the summation. 
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Thus, in the present illustration, we can write 


35 
(15 95 41 

41 95 

2135 


and such a product form will always be possible. 


6. Comparison with the Elliptically Loaded Wing 


From the general results already given, it is clearly possible to find the 
distribution of circulation K corresponding to any assigned downwash distribution 
w. A natural question to ask at this stage is: how many horseshoe vortices should 
be employed to give a distribution of vorticity leading to a total load correct to 
some given degree of approximation? This can clearly only be answered for some 
known case; the simplest and best known is, of course, the elliptically loaded wing. 


If the circulation at the centre of an elliptically loaded wing of span 2s is K,, 
then® the induced velocity has the constant value K,/(4s), while the total load is 


L= 5 pVK,s. « 


Suppose we represent this case by a system of m horseshoe vortices, prescribing 
the downwash velocity w; of equation (1) (and (5)) as K,,/(4s) for each element. Then 


K, 


The total load generated by this system of vortices will be 


2s 


n 


pVK, 
1 


on use of (24) this becomes 


Ln=— (25) 


Accordingly, the ratio of the load derived from the horseshoe vortices to the true 
elliptic load is 


-1 


This quantity will now be evaluated. 
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Since pre-multiplication and post-multiplication by a row and column of units 
is equivalent to adding all the elements of the matrix, it is found from (19) and 


(26) that 
jot tat ret 


Re 


where 
+j-1). .. -(k-r+1) (27) 


i 
But = B,=B,+B,+....+8B;, and therefore 


B,=B,+(B, + B,)+(B,+B.+B,)+. 


=mB,+(m—1)B,+. Br, 


(m+1—r)B,. 
1 


r= 


The factor (m+1—r) appears in the denominator of the expression for B,; hence 


(= 
= (j-1)!(m-}j)! 
Thus 


since the sum of the terms involving r is unity by Theorem B (Section 4). 


finally 
and if the factor (k + j— m-— 1) is written as the difference between (kK - m) and (1 - /) 


this may be written 


Rn = (j- 1)!(m-j)! 
m? j=) 
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OOWN WASH 

SPECIFIED 

AT POINTS 


(b) TWO VORTICES (¢) THREE VORTICES 


A A A A A 


(4d) SIX VORTICES (e) TWELVE VORTICES 


FiGuReE 2. 


By Theorem C the first summation equals — mk and by Theorem D the second 
equals 4m(m—1)k(k+1). Hence 


R,,.= [ (am — 2k) mk ~mk (m - I) (k +1) ] 


m+l1 


on reduction. For this illustration — the symmetrically loaded wing — k will 
normally be taken as 4, so that then 


_m+i1 


= 2 
Rn (29) 


7. Some Comments 


The result (29) is of some interest: it shows, for example, that if ten horseshoe 
vortices are employed to simulate the continuous distribution, the load generated 
will be ten per cent too high. 


Figure 2 shows some comparisons of loading, with 1, 2, 3, 6, 12 horseshoe 
vortices to represent the elliptic loading (shown for convenience as circular), and 
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TOTAL LOAD EQUALS ELLIPTIC LOAD, 
DOWNWASH CORRECT AT POINTS A, 


k=O 1464 


A AA A A A A A 


FIGURE 3. 


with k=4 for symmetry. The correctness of the formula (29) is demonstrated in 
these diagrams. 


The result (28) can be used to find a value for k which gives R,,=1 for all 
values of m: it yields on solution for k in this case 


As would be expected, k becomes equal to 4 when m is indefinitely large. Some 
distributions with k as defined by (30) are shown in Fig. 3: they are for m=1 
(k=0-1464), m=3(k=4) and m=8(k=4). The very considerable asymmetry of 
the latter two cases is to be noted; it suggests that a choice of k other than 4 is 
unlikely to be satisfactory. To obtain the results of Fig. 3, it is necessary to 
evaluate the column {K;} as given by (24). On use of (19), equation (24) becomes 
(with n=m-+ 1) 


_ § (- 
(m+1-n'(r- pt\(k-r+j 


and on use of Theorems C and D the summation involving j may be evaluated, so 
that finally 


(m+1-r)!(r—1)! | 
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8 A Second Example: Anti-symmetric Loading 


The case of symmetric loading provides the most appropriate test of the 
representation of a continuous load by horseshoe vortices: it is, however, desirable 
to show that the order of error demonstrated by equation (29) is not peculiar to 
this case. Accordingly, the rolling moment due to anti-symmetrical load is now 
examined. 


Consider a wing of span 2s on which the circulation is defined by 
K= - K, sin 20 
with yv=-—scos@, 
In the usual way, the induced downwash is obtained as 


K, Kyy 
w= - —cosd= . 


The total rolling moment L» is 


Es= = VAS. 


4° 


The loading (32) will now be represented by a finite set of m horseshoe vortices 
as before, with the induced downwash prescribed at the centre of each vortex (that 
is, k=4) by (34). This involves the solution of equation (5) with 


K, 
(2j - 1—m), ; (35) 


that is, {Ki} = -m}. . . (36) 


| This expression provides the ordinates {K;/K,} for comparison with the continuous 
curve of (32); the results are plotted in Fig. 4 for the cases m=4, 5, 6. As regards 
the total rolling moment, this is clearly given by 


2s s 
S pVK. 
pVK; (2i—1- mm), 


that is, by 2pV i}. . (37) 


In view of (35) and (36) the ratio R,, of the rolling moment due to the m horse- 
shoe vortices to that derived from the continuous distribution may be written 
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ANTISYMMETRIC LOADING 


OOWNWASH CORRECT AT POINTS A. 


FIGURE 4. 


The evaluation of this result, with m general, requires some labour: only the 
general lines of the analysis are given here. Since U~' is determined by (19) with 
n=m-+1 and (38) yields 


j=l t=1r=1 


G- 


(r—1)(m+1-r) 


thus we may write 


m) jy! 


where. on cancellation of the factor (mm + 1 - r) appearing in (39), 
G+ 1- r) (r-1) 


2 
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Yo ERROR 


ANTISYMMETRIC 


SYMMETRIC 


ASYMPTOTE _ 


--------o 


— 


10 20 


FIGURE 5. 


The quantity D, is readily evaluated by use of Theorems B and C; we find 


D,=4 (j-4—4m). 


Thus finally, with a little rearrangement 


m m 


This summation, involving as it does (m+ 3) factors in the numerator, cannot 
be evaluated by means of the Theorems quoted; however, by the methods of Ref. 1 
it yields finally 


(m— 1) (m+ 1) (im + 2) 


R, 3 (40) 
m 


In Fig. 5 are plotted curves of R,, as given by (29) for the symmetric loading 
case and by (40) for the anti-symmetric case. Both curves show that there is con- 
siderable error unless a very large number of horseshoe vortices is employed. 
The anti-symmetric case is nearly twice as bad as the symmetric; this is probably 
due to the fact that each half-span has a continuous loading of one sign and approxi- 
mates to a representation of symmetrical load by 4m vortices. 


§. Distribution of Downwash due to Horseshoe Vortices 


It is apparent from the preceding investigations that the choice of the points 
at which the downwash is specified in a horseshoe vortex system can profoundly 
affect the resultant load. This suggests that the downwash due to a horseshoe 
vortex system varies considerably across the span; and it is clearly of interest 
to examine this question. 
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Consider an elliptically loaded wing, and let the span 2s be divided into m 
equal parts; if the part of the elliptic load carried by the @™ section is L;, the mean 
circulation for that section, K;, is given by 


L;=pVK; (2s/m). 


If the set of circulations K; determined in this way is now used to define a horse- 
shoe vortex system, the total load will be correct, and the downwash for the system 
will be given by equation (1), in which k is now regarded as an independent variable 
— effectively the spanwise variable. 


Values of K; have been found in this way for a number of values of m; the 
downwash distribution has been calculated from equation (1). Figure 6 shows 
the cases m=1, 2 (which are identical, since for m=2 the two central trailing 
vortices cancel) and m=8. The downwash is expressed as w/w,, where w, is the 
constant value for the continuous elliptic load. 


As each (finite) trailing vortex is approached, the downwash tends to infinity; 
it is apparent that there is a very large difference from the constant value (w/w,=1, 
shown dotted) appropriate to the continuous load. Near the centre of the wing the 
downwash is roughly constant but, as the tip is approached, the variation becomes 
exceedingly rapid; presumably this is due to the increasing strengths of the trailing 
vortices toward the tip. 


In view of Fig. 6, it is perhaps remarkable that an arbitrary choice of position 
for specification of the downwash, such as k=4, should give a load of anything like 
the right order, either in magnitude or distribution. However, it is worth noting 
that in the case m=8, the intersections of the downwash lines with w/w,=1 occur 
very near the centres of the bound vortices, except for the tip vortex, where the 
intersection is noticeably nearer the wing tip. 


10. Practical Applications 


In the foregoing analysis comparisons have been made between the loads due 
to continuously varying distributions of vorticity and those due to horseshoe vortex 
systems; they have been related solely by the induced downwash. In practice, of 
course, the loads are usually only partly determined by the induced downwash, and 
our earlier results are only partially applicable in consequence. A brief examina- 
tion of this question is now attempted. 


Let us consider again the elliptically loaded wing, of elliptic plan form, un- 
twisted, and at geometric incidence z. The local load per unit length then gives, 
in the usual notation, 
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where K,, c, relate to the centre of the wing. 


But w=K,/(4s); 


ac, K, 
hence 4Vaa= Bs +1) 
and 
which, with a= 2r, yields the well-known result | probl 
2 
\Vcaz= 1)K. ‘ (22 
it ma 
The coefficient of K consists of two terms; the first, 2/A, derives entirely from 
the induced downwash, and the second from the geometric incidence. Equation 
(42) is usually taken to mean that the lift curve slope a is reduced by the factor 
A/(A +2); in reality, the geometric incidence is reduced by downwash from the 
value z to the effective value 
since 
Ww 2 Ais 
so that the induced downwash angle is 2/(A + 2) times the geometric incidence. 
The representation of the continuous elliptic load by a finite set of m horse- 
shoe vortices is now considered. The equation corresponding to (41) is 
hors 
with w; determined by (3). so that . 2/ : 
ist 
w= 
Accordingly, the equation corresponding to (42) becomes 
_fi = ] 45) Eq 
where c/c, is a diagonal matrix, with elements c;/c,. 
tert 


Now when A is indefinitely large, the first term on the right vanishes, and K 
is then determined wholly by the geometric incidence; the total load due to the 
horseshoe vortex system is then the same as that for the continuous elliptic load. 
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TABLE I 


ELLIPTIC WING, m=6: VALUES OF Rm (A). 


Equation (45) 1-167 1-102 1:034 1:0 
Equation (47) 1-167 1-105 1-037 1:0 


' On the other hand, when A becomes very small, the first term on the right of (42) 
/ and (45) becomes dominant. The evaluation of K then becomes precisely the 
problem dealt with in the earlier part of this paper, which led to the result (29). 


42) If (29) is regarded as implying that the elements of K are, to a first approxima- 
| tion, each greater by the factor R,, than those corresponding to the elliptic loading, 
it may be inferred that (45) is approximately equivalent to 


om 
ion 
tor 
since when A is large this gives the same value for K as the elliptic load, while when 
13 A is vanishingly small K is increased by the factor R,,. For A finite and non-zero, 
_ the new ratio R,,(A) is then, on comparison with (42), 
ARn A+2 
2 
4) 
When R,,~1 this shows that the fractional increase R,,—1, calculated for the 
horseshoe system when it is wholly defined by downwash, must be factored by 
2/(A +2) approximately to give the correct fractional increase R,,(A)—1; the factor 
is the proportion of geometric incidence due to downwash, as shown by (43). 
The formula (47) has been tested by calculation for the case m=6, for which 
R,=7/6. 
) — Equation (45) was solved for the circulations K;, from which the load was found; 
_ Table I compares the results obtained by use of the approximate formula (47). 
The formula (47), when equation (29) is used to define R,,, gives explicitly in 
terms of A and m 


(A+2)(m+])_ 
RAA)= (A +2)(m+1)-2° 
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Conclusions 


If a set of horseshoe vortices is to be used to represent a lifting surface, the 


foregoing investigations suggest the following conclusions. 


(a) 


(b) 


(e) 


If the circulations are to be wholly determined by the downwash induced, it 
is best specified at the centre of each vortex, except possibly near the tip. 


With the hypothesis of (a), a considerable error will result unless a large 
number of vortices is used: if the error is not to exceed 5 per cent, then not 
less than 20 vortices are required for a symmetric distribution, or 40 vortices 
for an anti-symmetric distribution. 


Correction factors based on equations (29) and (40) could be applied. 


If the circulation is determined by a geometric incidence modified by down- 
wash, then aspect ratio will play a part and the correction to be applied is 
given by equations (47) and (48); if the error is to be not greater than 5 per 
cent, the number of horseshoe vortices should be greater than 


With the hypothesis of (d), correction factors based on (48) could also be 
applied. 
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A Random Distribution of Gusts 
Corresponding to a Measured Frequency 
of Vertical Gusts 


M. FINE 
(Rotol Ltd.) 


SumMmMaRY: The paper describes a method of estimating the frequency of gust 
components in the plane of a propeller from measurements of vertical gust 
components, which are those generally recorded. By this means additional 
flight recording is avoided and the danger of under-estimating the frequency 
of gusts affecting the propeller is reduced. Random distribution of gusts in 
space is assumed. It is found that gust components in the plane of the propeller 
of 10 ft./sec. or more occur 2-51 times as frequently as vertical components 
of 10 ft./sec. or more. For speeds greater than 10 ft./sec., the ratio is higher. 


Introduction 


Accelerometers in cruising aeroplanes generally record vertical components 
of gusts. The propeller, however, is stressed by all components in its plane, rather 
than by vertical ones alone. Lateral gusts will not be recorded by an accelerometer 
measuring vertical components; hence available records of vertical components 
under-estimate the frequency of gusts affecting the propeller. The problem is to 
estimate the incidence of such gusts from the measured vertical components, thus 
avoiding additional flight recording. 


NOTATION 


The term “vertical gusts” is used here for “vertical components of a gust”, and 
“propeller gust” for the “component of a gust in the plane of the propeller.” 


u vertical component of gust velocity (ft./sec.) 
component of gust velocity in propeller plane (ft./sec.) 
inclination of direction of v to horizontal (radians) 


number of vertical components with velocity u or more 
for every 1,000 of 10 ft./sec. or more 


number of propeller gusts with velocity v or more for 
every 1,000 of 10 ft./sec. or more 


_ frequency of propeller gusts of 10 ft./sec. or more 


~ frequency of vertical gusts of 10 ft./sec. or more 
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number of vertical gusts recorded with velocity 1 
or more if there are 1,000 d6 propeller gusts of 
10 ft./sec. or more in the directions between 6 and 
6+d6 


A,B,2z,8 constants used in the expression for g 
G(v) function of v7 used in the expression for g 
p,r,w,x,y,& auxiliary variables 
h,H auxiliary functions 
a, auxiliary parameters 
c=(2/ x)". 


2. Assumptions 


The vertical gusts are described by Bullen; 
F (u)=27,800 exp (—0-34411u) + 878-2 exp (—0-20816u) . » 


fits measured results well. (Here “vertical” means both “up” and “down.”) Gusts 
are assumed to be randomly distributed in space and the distribution which gives 
vertical components agreeing with relation (1) is found. 


3. Mathematical Formulation of the Problem 


A gust with a component normal to the propeller is resolved into this com- 
ponent and one in the plane of the propeller. The components normal to the 
propeller are neglected for the time being; clearly the remaining components are 
distributed randomly. The following argument is applicable to any random 
distribution of vectors in a plane with records of components in one direction. 


Figure |, in the plane of the propeller, shows a component of velocity v inclined 
at 4 to the horizontal. The vertical component is v sin @=u, i.e. v=u/sin 6. 


Consider an element defined by dé and propeller gusts in the corresponding 
range of direction, i.e. these propeller gusts have no components at right angles to 
the direction 6. If there are 1,000 dé propeller gusts with speed of 10 ft./sec. or 
more in the range given by dé then there are g(v)d@ propeller gusts with speed 
v ft./sec. or more in this range, and g(u/sin 4) d@ propeller gusts with vertical 
components of u ft./sec. or more. 
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The number of vertical gusts of u ft./sec. or more recorded by the accelero- 
meter is 


«/2 


N (u)=2 g(u/sin 6) dé. . (2) 


0 


The numerical factor is 2 because F, in relation (1), includes “up” and “down” gusts. 
The corresponding number for vertical gusts of 10 ft./sec. or more is 


N (10)=2 | g (10/sin 4) dé. ‘ & 
For every 1,000 recorded vertical gusts of 10 ft./sec. or more there are 
1,000 N (u)/N (10) vertical gusts of u ft./sec. or more, and so 
2/2 


(2.000/N (10)} | (u/sin 6) (u) 


0 


= 27,800 exp (— 0:34411u) + 878-2 exp (—0-20816u). 


(4) 


It is necessary to find g from (4), with N (10) given by (3). g is the function for 
propeller gusts corresponding to F for vertical gusts. If g is known, R can be 
found. 


With 1,000 dé propeller gusts of 10 ft./sec. or more corresponding to each 
element d@, the total number of such gusts is 


2/2 


2 | 1,000d6 = 1,000. 


0 


The accelerometer has recorded N (10) vertical gusts of 10 ft./sec. or more. 


Hence R= 1,0007/N (10). 


The method of solving equation (4) is outlined in the Appendix. 


4. Results 
The solution of equation (4) is 


AG (av) + BG (fv) _ 
AG (102)+ BG (108) ° 
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g (v)=10° 


u 
d 
2/2 


where A =27.800, 
B=878-2, B=0-20816 } 


and G (x)=(2/7) xK, (x), . (9) 
K, (x) being a Bessel function* of imaginary argument®: *: ”; 
R=(z/2,000) [AG (10z) + BG (108)]. ‘ (10) 


Figure 2 shows the variation of g with v and the variation of F with u. The value 
of R is 2:51, so that recordings of vertical gusts fall well below the number of 
propeller gusts. 


Components normal to the propeller have been neglected; by symmetry they 
occur as often as vertical gusts. 


F (u) MEASURED, VERTICAL, u 
‘ 


| 
iN \ | 
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| | | | | 
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| | | | 
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Appendix 


METHOD OF SOLVING EQUATION (4) 


n/2 
If G (v) satisfies G (v/sin 6) dé=exp v), 
0 


then the solution of equation (4) is as stated in equations (7), (8) and (9). 


oo 
Write G(v)=c h(p) cos pudp, 
0 


where 
by Fourier’s integral theorem 
oo 
h(p)=c | G (v) cos pudv. 


0 


Substituting from (12) in (11), putting 


p=r sin 6, 


and inverting the order of integration gives 


oo 
c | H (r) cos rudr=exp (- v), 
0 


where H (r)= f h(r sin @) sin 6d6. 
0 


Applying Fourier’s integral theorem to (16) gives, with (17), 


co 
H (r)=c | exp (- v) cos rudv=c/(1+r*)= h(r sin sin 6d6 
0 


0 


Expanding (1 +r?)-! and A (rsin 8) in powers of r? and equating coefficients gives, 
with A(rsin #)=a,+a,r* sin? @+a,r' sint 6+ .. 


a,=c, a, = 3c/2, a, = (3 x S)e/(2 x 4), 


and h(p)=c (14+ p?)"?%; 
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(11) 


(12) 


(13) 


(14) 


(15) 


(16) 


(17) 


(18) 


(19) 


(20) 


(21) 


| 
) 
f 
= 


whence the step from (15) to (16) and the result (21) may be justified as follows. 


exp (-v)=c [| h (p) cos (pv /sin 0) dp | a8 
0 


{ h(r sin 6) cos rvdr sin 6 dé. 


0 


x, 
[ h (r sin cos rvdr sin dé = [ h(r sin @) sin | cos rvdr. 
0 xX, x, 


With h given by (21), the last integral in square brackets is (1 +r*)~! if r is complex 
and |r|<1. In the strip | .Y (r)| <1 the integral is an analytic function of r and, by 
continuation, is (1 + r?)-?. 


Hence the last double integral is 


x, x, 


| (1 cos rudr r-*dr (in absolute value), 


1 


which is smaller than any assigned quantity for sufficiently large X,, X,, and the result 
follows. 


From relations (12), (13) and (21), 


0 


Equivalent integrals were discussed by Glaisher®) in 1881, but his work was not 
seen by the author until after the solution had been found independently’. Using 
Glaisher’s results gives the solution more quickly than otherwise; Ref. 6 uses rather 
less advanced analysis, yet provides some interesting results. 


Glaisher shows that, with a slight change of notation in his work, 


=x” cos agdé 
(x? 


0 
is a solution of the Riccati equation 


m(m-I)y 


x? 
if y is finite. 
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Putting €=xp and ax=v in (23) gives 


y= cos pudp 


0 
Hence if m=1'5, 


and, from (24), y is a solution of 


m(m-—1)y 
= 


2 1 0-5 2 


This has as solutions Bessel functions of order y= m — 0-5, with imaginary argument iv. 
As v tends to infinity x tends to infinity and, from (23) and (28), both y and w tend 
to zero. This disposes of one of the solutions of (27). The behaviour of y and w 
when v tends to 0 gives the multiplicative constant and, for m=1-5, it is found that 


G (v)=(2/) vK, (v), ‘ (30) 


where K, is as defined in Ref. 2. 


This method. may be applied to find g(v) corresponding to other forms of F (u) 
and, in addition, to find the three-dimensional distribution of gusts. 
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The Influence of the Wake on the Flutter 
and Vibration of Rotor Blades” 


J. P. JONES, Ph.D., B.Sc.(Eng.) 


(University of Southampton) 


SuMMARY: An approximate theory is described which takes into account the 
influence of the wake on the aerodynamic derivatives of an oscillating rotor 
blade. It is shown that the most important parameter is the number of cycles 
of the oscillation which occur during one revolution of the rotor. For integral 
values of this parameter the flexural damping is shown to become very small 
at small angles of incidence. The flutter characteristics of a simple rotor blade 
are Calculated on the basis of this theory, and it is found that the wake causes 
radical changes in the critical speeds. The theory is in good qualitative agree- 
ment with experiment. 


1. Introduction 
1.1. Rotor BLADE BENDING FREQUENCIES 


Rotor blades derive much of their bending stiffness from the centrifugal force; 
in fact in the fundamental mode of a hinged blade the centrifugal force supplies all 
the stiffness. It follows that the natural bending frequencies increase with 
rotational speed, the increase depending upon the mode and the root conditions. 


Experimentally, it is found that 


where z is a constant whose value depends upon the mode and the root conditions. 


Approximate values of z for hinged and encastré blades are: — 


Hinged Blades Encastré Blades 


Fundamental 1-0 Fundamental 
Ist Overtone 6:0 Ist Overtone 
2nd Overtone 12-0 


*The material in this paper is based on a report submitted to the Department of Scientific (reve 
and Industrial Research in June 1954, which was later extended and circulated as A.R.C. Report It is 
18,173, in January 1956, The present paper was received for The Aeronautical Quarterly in strai 
October 1957. 
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NNATURAL FREQUENCY 


/ 
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REV/SEC. 


FiGureE |. Variation of first bending frequency with rotational speed for a hinged blade. 


A typical curve of natural frequency (cycles/sec.) against rotational speed 
(revolutions /sec.) for the first overtone mode of a hinged blade is shown in Fig. 1. 
It is assumed that f, =8-65 c./s. and z=6-0. Superimposed on the diagram are the 
Straight lines “3,” “4,” “5,” etc. The enclosed number denotes the slope of the 
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VIBRATORY 
STRESS 


(crctes /sec) 
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9 


NATURAL FREQUENCY 


2 


RESONANCE BENOING MOMENT 
STEADY-STATE BENDING MOMENT 
BENDING FREQUENCY 


° 


70 


ROTATIONAL SPEED (R.RM.) 


FiGURE 2. Hovering case. Hinged blade. Blade angle 12°. Variation of stress with rotational 
speed"), 


line in cycles per revolution. Thus, at the rotational speeds corresponding to 
the intersection of any of these lines with the frequency curves, the natural frequency 
will be such that an integral number of cycles are completed in one revolution. 


1.2. EXPERIMENTAL RESULTS 


Hirsch®’ and Meyer’? have recently published the results of some experimental 
work on the stresses in rotor blades, and both found unexpectedly large fluctuating 
bending stresses, with frequencies varying from one to ten cycles per revolution. A 
diagram of the type shown in Fig. 1 was constructed for each blade, and it was found 
that the large stresses occurred at rotational speeds such that one or more of the 
natural bending frequencies was given by 


frequency = integer x rotational speed. 


On the basis of the usual aerodynamic theory, these vibrations cannot be 
explained. The results certainly suggest that the blades were in resonance with 
exciting forces having frequencies of one per rev., two per rev., and so on, but if so, 
what was the source of excitation? In forward flight, the lift varies periodically as 
the blades rotate, but theory predicts only very small higher harmonic components, 
which reach a maximum at large », whereas the measured stresses were largest 
when « was small, i.e. in the transition region.* 


Daughaday and Kline“ have put forward the theory, based on a paper by 
Mangler and Squire, that the excitation is produced by the variation of induced 


*The “tip-speed ratio” » is defined as V,/QR, where V,, is the forward speed of the helicopter 
and R is the tip speed due to rotation. 
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velocity around the disc. By induced velocity is meant the downwash at the rotor 
disc due to the trailing vortex system. The agreement between theory and experi- 
ment was not good, due possibly to the difficulties of obtaining a reliable estimate 
of the downwash, especially near the blade tips, but more probably to the analysis 
being incomplete. In forward flight the induced velocity is certainly not constant 
around the disc, and a forced oscillation is possible; but it appears from Meyer’s 
results that the vibrations may also be experienced on a hovering rotor (in which 
case the induced velocity is constant around the disc), and in the full scale tests of 
Ref. 2 the higher harmonic vibration was first noticed during a ground run. To 
investigate this phenomenon Meyer carried out a comprehensive series of model 
tests, on hovering rotors, using both hinged and encastré blades. Figure 2 is 
reproduced from Meyer’s report and shows (a) the variation of natural bending 
frequency with rotational speed, and (b) the variation of the ratio (resonance b.m.)/ 
(steady state b.m.) with rotational speed. The most important features of the 
results are now described. 


(i) Within the range of the experiments, the occurrence of the vibration was 
independent of the mode of oscillation and of the absolute value of the frequency. 
The fundamental (encastré) and first overtone (hinged) modes were excited at 
frequencies up to eight cycles per revolution. The peak vibratory stresses were 
quite large, being as much as 50 per cent? of the steady bending stress in some cases. 


(ii) It was found that the amplitude of the vibration increased with increase 
of incidence but, in particular, the oscillations occurred at zero incidence with blades 
of symmetrical section, i.e. when there was no thrust, and therefore no trailing vortex 
system. The amplitude did not remain constant, but grew and decayed in an 
irregular manner. 


(iii) Since a wide range of blade stiffnesses and frequencies was used, the 
possibility of the supporting structure playing an important part in the vibration 
could be eliminated. 


(iv) It will be seen from Fig. 2 that the vibratory stress increases very rapidly 
as the frequency ratio v (number of cycles per revolution) approaches an integral 
value, but that it never becomes zero. Therefore it must be assumed that a small 
periodic force was present at all speeds, produced either by fluctuations of torque 
or by slight eccentricity of the shaft and bearings. Adjustments to the balance and 
alignment of the rotor produced only slight changes in the amplitude of vibration. 


The stress-speed curves are similar to those obtained if a lightly damped 
structure is in resonance with a periodic exciting force, but simple aerodynamic 
theory (i.e. assuming that the oscillation does not affect the lift-curve slope) predicts 
heavy damping of flexural vibrations. The possibility that the exciting force itself 
increases rapidly at certain rotational speeds may be ruled out, because with differ- 
ent blades the peak amplitudes were obtained at different rotational speeds, for the 
same frequency ratio. 


{The vibratory stresses in Fig. 2 are the peak-to-peak values, i.c. twice the amplitude. 


August 1958 261 


Thus it was concluded that the vibrations were either (a) self-excited, or (b) 
forced, but simple aerodynamic theory was invalid. A brief review of the evidence 
shows that if (a) is the explanation, then the vibrations were so much divorced from 
the usual types of instability that they must be regarded as completely new 
phenomena.* Whirling would involve some form of coupled motion between the 
rotor shaft and the blades, with perhaps as many as five critical speeds quite close 
together, but with the same mode of blade bending at all five speeds. For similar 
reasons the usual types (i.e. fixed-wing) of classical flexure-torsion flutter can be 
eliminated. The vibrations occurred in one mode only, at the natural frequency of 
that mode, and the several critical speeds were quite close together. The construc- 
tion of the blades was such that both the flexural and inertia axes were substantially 
on the quarter-chord line; hence the aerodynamic and inertia couplings were very 
small. It is true that the rotation introduces inertia and velocity couplings which 
are absent from fixed wings, but in general these are small, especially at small 
angles of incidence. Stall flutter may be ruled out at low angles of incidence (but 
see Section 2.2). 


The alternative (5) is the more likely explanation, i.e. that simple aerodynamic 
theory is invalid. This implies that for integral values of v the damping becomes 
very small. It is known from oscillatory aerofoil theory that the flexural damping 
is less than that obtained if the effect of the wake is neglected. But the frequency 
parameter over the outer sections of a blade—where most of the lift is developed— 
is small, and the overall reduction of the flexural damping would not be more than 
20 per cent, which is by no means sufficient to explain the sharp increase of 
amplitude. There is, however, a significant aerodynamic difference between fixed 
and rotating wings, namely the form of the wake. Behind a fixed wing it is 
essentially a plane surface, extending to great distances, so that the remote parts 
have little influence on the flow past the wing; but the wake of a rotor blade or 
airscrew is a helical surface and the blade actually passes over its own wake. 


It is now suggested that the wake is primarily responsible for the hovering 
vibrations and, in the remainder of the paper, an approximate theory is developed 
from this idea. 


NOTATION 
f/2= natural frequency of rotating blade in a given mode 


f,/2= natural frequency of stationary blade in the given mode 
angular velocity of the rotor 


t time 
frequency 


*The possibility of the vibrations being essentially non-linear cannot be ruled out, because the 
amplitude did not grow without limit, and the wave form was not periodic. But if the blades 
had slightly different natural frequencies, then the wave form, and the fact that the amplitude 
had a maximum value at a certain speed, can be explained equally well by a forced oscillation. 
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(6) v=p/Q) frequency ratio (cycles per revolution) 
nce N integral values of the frequency ratio 
pe 6 fractional values of the frequency ratio 
i ¢ vorticity shed into wake at time ¢ 
Ose 2g vorticity shed into wake at earlier instant 
ilar 6 angle measured from trailing edge 
| 7 v=Or 
uc: c chord 
illy r radius of blade element 
ery 
6m element of mass 
all tip radius 
Dut r=c/R 
co distance of inertia axis aft of the quarter-chord line 
nic om» minimum value of o to produce flutter 
“ ck radius of gyration of blade about the spanwise axis 
cy ck, radius of gyration of mass-balanced blade about the spanwise 
axis 
an B blade flapping displacement (downward flapping positive) 
of x blade pitching displacement (nose-up pitching positive) 
. ly moment of inertia of the blade about the flapping hinge 
rts I, blade pitching moment of inertia 
of Ky root spring stiffness 
,; flutter speed coefficient 
p density 


n=2npcR* /(8/;). blade inertia number 


x= —4c cos %, co-ordinate of a point in the wake measured from 


mid-chord 
x, co-ordinate of a point on the aerofoil measured from mid- 
chord 
= 4, 
n=€-€, 


V=Or, free stream velocity at representative radius 

w’ reduced frequency at representative radius (=pc/2r) 
. C (iw’)=F +iG, Theodorsen’s Function 

X, (iw’)=C (iw’) I, (iw’) + 1—C (iw’) 1, (io )=X +i¥ 
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J, (o’), J, Bessel Functions 


(iw’), 1, (iw’) modified Bessel Functions 
I,’ =2 [cosec » - (1 —C) 4 cot 9] 
I,” =2 [4 cot cosec + iw’ sin 
[= -2sing+4 cot ¢+4iv’ sin 2¢ + iw’ sin » 
= —2 sin no +io’ [sin {(n+ 1) /(n+ 1)—sin {(n—1) /(n- 1)] 
K,”=c sin 
K,=4c sin ¢ +4 sin 29 
K,,= 4c [sin {((n+ 1) /(n+ 1) - sin {((n— 1) 9} /(n- 1) 
W’,. downwash due to vortex sheets in wake below aerofoil 
P,Q, R,S,X’, Y’ functions of ’, h and 4 defined in the text 


129 "22° "32 \ derivative coefficients 


Myo, Moo, 


2. Proposed Explanation 


2.1. EFFECT OF THE WAKE ON AN OSCILLATING BLADE 


One of the many difficulties of rotor aerodynamic theory lies in the fact that 
there is usually more than one blade, and the mutual interference of their wakes and 
“bound” vorticities cannot be ignored. In vibration theory, a further complication 
arises from the elastic and inertia couplings which exist between the blades, through 
the root and rotor shaft. Thus, in order to keep the theory as simple as possible, 
only a single blade rotor is considered, but it is recognised that this may introduce 
considerable errors. 


The blade is assumed to be rotating with constant angular velocity 2 and to be 
oscillating at a frequency of p radians per second. The oscillation is assumed to 
have existed for sufficient time for the wake to extend to infinity below the rotor. 
and for the amplitude to have become steady. Then the vorticity which is shed 
from the trailing edge will move along a helical path relative to the blade and, after 
one or more complete revolutions of the rotor, will be directly beneath the trailing 
edge. Since the inflow velocity through the disc is usually very small compared 
with the tip speed of the blade, the pitch of the helix is small, and the vortex sheet 
beneath the blade will be quite close to it. Thus the blade is disturbed by vorticity 
shed in previous revolutions, and if these disturbances were all additive they would 
tend to sustain the vibration, and the effective damping would be reduced. The 
induced effects of the wake will be a maximum if the vorticity below the blade is 
the same after each revolution, i.e. provided that an integral number of cycles of the 
oscillation are completed per revolution. 
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There are three main objections to the preceding theory : — 


(i) The vorticity does not remain absolutely stationary in space, but tends to 
“slide” along the surface of the helix. Thus, after each complete revolution the 


vortices of equal sirength would be shifted relative to each other, and the maximum 


induced effects would be obtained at a slightly different frequency. Meyer’s results 
show that this is not critical, since the vibrations occurred over a range of rotor 
speeds on either side of that which produced maximum amplitude. 


(ii) The vorticity diffuses with time. This is a more serious objection which 
obviously depends upon the rate of diffusion, but those parts of the wake in which 
an appreciable reduction of strength has occurred will be several chords length 
below the blade; hence the total error may not be large. 


(iii) The influence of the wake must decrease with increase of incidence, but 
the experiments showed that the amplitude of vibration increased with incidence. 
It is not possible at present to give a precise reason for this, but there are two 
probable explanations : — 


(a) If the vibrations were excited by small variations of the driving torque, 
then the amplitude of the forcing function would increase with incidence. 
The torque fluctuations cause small changes in the rotational speed, and 
hence fluctuations of the total lift, but the steady lift coefficient increases 
with incidence; hence the oscillatory lift component due to changes in the 
rotational speed will increase with incidence. 


(b) At the greater of the two incidences (0° and 12°) investigated, large areas 
of the blade may have been stalled, especially since the Reynolds number 
was very low. There is some confirmation of this from the Cornell 
Aeronautical Laboratory (C.A.L.) experiments (see Section 2.2), which 
were conducted at a greater Reynolds number, but stall flutter was 
obtained at a blade angle of 12° (in some cases 9°). 


2.2. EFFECT OF THE WAKE ON FLUTTER 


A further problem which arises from the foregoing theory concerns the effect 
of the wake on the flutter derivatives of a rotor blade. Assuming that the flexural 
damping can become very small at certain frequency ratios, there is no a priori 
reason why the remaining derivatives should not change radically at some value of 
v, and a recent paper from the Cornell Aeronautical Laboratory suggests that this 
may actually happen. At about the time this theory was being worked out, further 
experimental work was done at the C.A.L. on the flutter of a single blade model 
rotor and, through the courtesy of Mr. H. Daughaday, a copy of the report has been 
made available to the author. 


The model blade was free to flap about a hinge on the axis of rotation, and to 
pitch about a spanwise axis. The pitching motion was restrained by springs at the 
root, but in flapping the restoring moment was provided mainly by the centrifugal 
force. At the root an elastic coupling was introduced between flapping and pitch, 
and the pitching moment of inertia and chordwise c.g. position could be varied by 
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means of internal balance weights. Three blades were used, one of very. high 
torsional and bending stiffnesses, so that the motion was substantially that of a 
rigid blade, and two others of equal (low) bending, but different torsional stiffnesses. 


Two principal types of flutter were recorded, namely, classical flapping-pitch 
and a much milder vibration, which occurred only at small angles of incidence and 
at a frequency ratio of approximately 2:0. To account for the latter type of flutter, 
Daughaday and Kline put forward an explanation similar to that already given, but 
the modes of vibration, and the conditions under which it was obtained, differed 
considerably from those of Ref. 3. No single-degree-of-freedom vibrations were 
found, it was always necessary to provide some coupling (inertia or elastic) between 
flapping and pitch, and the frequency ratio corresponded to that of the coupled 
motion (v ~ 2:0). 


In some cases it was found that a small change in root spring stiffness caused a 
jump in the frequency ratio, from v=1:8 to v=2:2 and even from v=1°7 to 
v=2°6, i.e. for a small change in one of the parameters it was possible to jump 
from one type of flutter to another. 


It is possible, however, that the method of recording used by Meyer was more 
sensitive—or at least more revealing—than the C.A.L. technique, since the latter 
merely measured the rotation of the blade about its hinges, whereas Meyer’s models 
were extensively strain-gauged along the span. On the other hand, it is likely that 
torque fluctuations were less on the C.A.L. models, since the rotor was driven 
through a gear box, flywheel and a long shaft of appreciable moment of inertia. 


The elasticity of the blades was shown to have an important influence on the 
critical speeds, but it appeared that, for the configuration chosen, flexure was of 
more importance than torsion. (But see Section 5.4). 


In several of their investigations into classical flutter, Daughaday and Kline 
obtained the curious result that v remained virtually constant at about 1-50, inde- 
pendent of the incidence or coupling. This was not in agreement with calculated 
solutions, based on quasi-steady theory (i.e. the influence of cast-off vortices is 
completely neglected), which showed a substantial variation of v, and it must be 
concluded that the aerodynamic theory was at fault. 


In the next section an approximate theory is described which allows for the 
effect of the wake on the aerodynamic forces, and which gives results in good 
agreement with experiment. 


3. Theory 
3.1. AN APPROXIMATE REPRESENTATION OF THE BLADE AND ITS WAKE 


The forces and moments on an oscillating rotor blade may be calculated by a 
straightforward extension of the theory of oscillating wings of finite span, but this 
would be a very laborious process and some form of automatic computer would be 
essential. The greatest difficulty lies in the evaluation of the downwash integrals 
for a helical wake. Again, until it can be established that the wake does have an 
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important influence on the derivatives, and that the assumptions involved are 
acceptable, the labour necessary to obtain an exact solution cannot be justified. 
However, it is possible to cut down the labour considerably if the three-dimensional 
problem can be reduced to one in two dimensions, but which preserves the essential 
features of the former. 


This is a common process in aerodynamics, but especially in propeller theory, 
where it has been used by Prandtl to estimate tip losses, and by Joukowski to study 
the “cascade” interference between blades. The helicopter rotor is especially 
suited to this type of reduction, since the aspect ratio is high and the blade may be 
regarded as an oscillating aerofoil of infinite span. 


Now the radius of curvature of the outer sections of a blade—where the 
relative velocity is greatest—is large; hence, for a considerable angular rotation the 
outer blade elements may be assumed to move in straight lines. The linear distance 
travelled will be large compared with the chord, and since the pitch of the helix is 
very small—only a chord’s length or so at moderate incidence—it can be assumed 
that the wake is plane for several chords’ length ahead of, and behind, the aerofoil. 
Thus the successive turns of the wake which are beneath the blade can be regarded 
as equally-spaced, plane vortex sheets, parallel to that attached to the trailing edge 
of the blade. The sheets are not horizontal, but are inclined at the mean pitch angle 
of the helix, i.e. parallel to the relative wind. But if the aerodynamic forces are 
calculated on the assumption that the blade is vibrating about zero mean incidence, 
then the simplified configuration shown in Fig. 3 is valid. The relative wind 
(velocity V) is assumed to be horizontal. and the vortex sheets are spaced a distance 
hc/2 apart, where c is the chord and / is a function of the mean incidence of the 
non-vibrating blade. It is convenient in the analysis to assume that the vortex 
Sheets extend to infinity in both directions; hence it will be seen that the problem 
has been reduced to that of finding the forces on an oscillating, two-dimensional 
aerofoil, but with additional vortex sheets extending to infinity ahead of, behind, 
and below it. 


So far the distribution of vorticity in these sheets has not been considered, but 
if v is an integer, then the free vorticity at a distance x from the aerofoil (measured 
positive to the rear and from the mid-chord) will be the same in each sheet, but if v 
is not an integer it is necessary to consider the vorticity distribution in the 
helical wake. 


3.2. THE DISTRIBUTION OF VORTICITY IN THE WAKE 


Assume a system of cylindrical polar co-ordinates (r, 6, z) to be associated with 
the rotor; z is measured parallel to the axis of rotation and @ is measured from the 
position of the blade trailing edge at time ¢. If the oscillations are simple 
harmonic, then the vorticity shed into the wake at time ¢ and radius r is 


where « (r) is the amplitude of the vorticity at radius r. 
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Ficure 3. Aerofoil and vortex sheets. 
Note: O is not the origin for measurement of x. The mid-chord is the origin in Section 4.1. 


Then the vorticity cast off at this radius, but at some earlier time, is 


If 6=2nr, n=1, 2, 3 , then this will be the magnitude of the free 
vorticity in successive sheets directly below the trailing edge. 


Also, putting p/2:=v=N+4, where N=1, 2,3... ., and 0<6< 1-0, then 


i (N+8) 2nz 


=e@e"e 


i.e. in neighbouring sheets the points of equal strength are displaced relative to 
each other by an amount which depends upon the frequency ratio. 


This result can be transferred directly to the approximate representation if it is 
assumed that the curved path of the blade can be represented by a straight line. 


Put @=(x -c/2)/r=c(€-1)/(2r). 
Then 


=e, (re 


where ¢,(r) is the peak value of the vorticity at the trailing edge. 
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For the sheet which is attached to the trailing edge 1 < £ < 00, and for the 
remaining sheets - 00 < £ < 0%, 

If the blade path is straight, then the strength of the vorticity at any point in 
the sheet attached to the trailing edge (n=0) is 


and the free vorticity at the same value of € in the n™ sheet is 


The free vorticity at the trailing edge can be obtained from the usual theory of 
oscillating aerofoils. 


4. Theoretical Analysis 
4.1. FORMULAE FOR THE LIFT AND MOMENT COEFFICIENTS 


This section is based on the theory of oscillating aerofoils as described by 
W. P. Jones. 


If z, is the downward displacement of any point x, on the chord at time f, it is 
necessary to satisfy the boundary condition of tangential flow at the surface, i.e. 
the equation 

dz 
=w(x,,Z,, 0) . : (9) 
Ox, 
where w(x,, z,, f) is the downward induced velocity due to the bound and free 
vorticity. 


If the frequency of the oscillation is p/(27) c./s., then we may put 


i 


and we assume that 
the bound vorticity y= Te” 
free vorticity == Ee” 


ipt 


total circulation k= Ke 


and the downwash w=W e”™ 


[ is written in the form 


P=V(C,@ 
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where the functions I",’, I',”, I’,, and so on, are defined in the Notation, and the 
constants C,, C,, and so on, are to be determined. 


It can be shown that for points on the aerofoil 


and for points in the wake 
+1 


K=Se « .« & 


Hence, from each of the functions I",’, I’,”, and so on, the corresponding functions 
K,’, K,”, and so on, can be derived (see Notation). With the exception of 
K,’, the K, are zero in the wake. 


The downwash w (x,, z,, t) at the aerofoil is obtained from 


= | KS 


dx \(x—x,)? + (nhe - J 
and z, —> 0. 


For an aerofoil of infinite span it is found that 


W,’=0, W,”=10, W,=4+cos¢, W.=cosn, . . (16) 


but in the present case W,’ 0 because of the induced effects of the additional 
sheets. Let the induced velocity of these sheets be W’,,; then 


+00 
1 


(18) 


X,=C 1, (io) + [1 (io’)] 1, (io’) 


iw’é 
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and, if we put €-,=n, 
+00 

j 


n=X 
1) 


n= {P+ 


In the Appendix it is shown that 


+00 
lw 7 


{n? + 


iw'nh 


dyn= — ite 


Hence W’,,=zw’X,e Se”. 


n=1 


where C= 0h + i276. 


n 


n=1 


ys TO i 
W’, 


where R= 
~ {[e** cos 276 — 1]? + sin 


—e*" sin 25 
{[e** cos 278 — 1}? + sin 


Now =e™ — (w’ cos +i sin (w’ cos 


and it may be shown that 


cos (w’ cos ,)=J, (w’)— 2 {J, (w’) cos 20, —J, (w’) cos 49, +... 


sin (w’ sin ¢,)=2 {J, (w’) cos ¢, —J, (w’) cos 39, +... }, 


where J, (w’), J, (’), and so on, are Bessel functions of the first kind. 


Because of the high aspect ratios of rotor blades, the reduced frequency is 
usually small (w’ <0:25). If attention is restricted to small ’, then, to a good 


approximation, 


(R+iS)=XJ, —2,Y’ cos, +i(YJ,+2XJ, cos¢,) (31) 


and 
,_ {X [e*" cos 25 — 1] + Y e* sin 278} 
where X’= cos 2nd — sin ‘ (32) 
271 
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y’= {Y cos 2x6 — 1] — X sin 275} zw’ 
cos 275 — 1]? + [e** sin 
Then at some point #, on the chord, the downwash is given by 
W (¢,,2,, 0=V [C, 1+ W’,,)+C, (4+c0s ¢,)+C, cos2¢,J+.... . (34) 


and the constants C,, C,, C., and so on, can be obtained from the relation (sce 
equation (9)) 


When the C,, are known, the downward force L and moment M, about the 
quarter-chord, per unit span, can be obtained from 


-L= 4pVe | sin do 
0 


pvc? sin {4 —cos do 
0 


M 


it is convenient to express C, and Cy in the form 


Putting Ci 


Cy +m,,x + +m,,c*x] V7}, (40) 
where z is the downward displacement of the quarter-chord axis and x is the 
rotation (positive nose-up) about the quarter-chord axis. 


The quantities /,., /,,, and so on, will be referred to as the derivative coefficients 
and in Section 4.2 they are evaluated for the case of a blade performing combined 
translational and rotational oscillations. The first number in the suffix refers to 
the type of derivative, i.e. displacement, velocity or acceleration, and the second 
number to the type of motion, i.e. translation or rotation. 


4.2. CALCULATION OF THE DERIVATIVE COEFFICIENTS 


Let the downward displacement z of the quarter-chord axis be 


— ipt 
4cz (41) 


and let 
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Then the tangential flow equation (35) can be written 
W/V =iw’ Z+ x, [iw’ (4--cos 1] ‘ (42) 
or C,(14+W’,)+C, G+cos 9)+. . . =iw’ Z+y, {iw’ . (43) 


Substituting for W’,. from (31), and equating the coefficients of the trigonometrical 
functions, gives 


C,=(P+iQ) [iw’ z+ x,(U+io’)] ‘ . (44) 
and C, = — iw’y, + 2J, (w’) [Y’P + OX’ —i (PX’ -- (45) 


and Q= (YJ, -X’J,}? (47) 


C,, C,, and so on, are zero to this approximation. 


Now, for a wing of infinite span, 
P= 1-0, Q=0, 


and, if the wing is performing translational oscillations alone, i.e. x,=9, 


Z | 
and C,=0. J 


C, is only different from zero if the wing is pitching, but in the present case 
when x,=0,C, 40. Thus the wake produces an apparent rotation of the aerofoil. 
Similarly, an apparent camber change is produced if more terms are retained in the 
series expansions for cos (w’ cos ,), sin (w’ cos @,), since then C,, C,, and so on, are 
not zero. 


Now from (36), (37) and (38) 


{ (iw’) + 


(49) 


Then, on substituting for C, and C, from (44) and (45), the following expressions 
are obtained for the derivative coefficients. 
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1,=FP-@QG - PGu'+ 2 ) (PX’ - 


4 
i+ (Y’P+QX’ )| 
- 20" (OF + PG) 


- Qo’ oJ, ) (PX’ - OY’) 


| 


M,,=0'Q +25, + (w’) (PX’ - OY’) 


m,,=4J, (w’) (PX’ - QY’) 


m,,=Qw’ +2J, (w’) QxX’ - (PX’ - oy’) | 


m,,=3| -P- 
(51) 


X’, Y’, P and Q are functions of incidence, reduced frequency and 6. No 
attempt has been made to provide comprehensive tables of these functions, but P 
and Q have been calculated assuming the values given in Table I. 


TABLE I 
h 
0°10 

0:0833 
0-250 
0:50 
0-75 
0:9167 
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| 43. DISCUSSION OF THE RESULTS 


In Table I, w’ was estimated on the assumption that the blade was moving in a 
| straight line with the velocity of an element at about 0-75 tip radius. The three 
values of A correspond to zero, small and moderate incidences, also at 0-75 tip 
| radius. No attempt has been made to deal with very large incidences. The 
j variation of the most important derivative coefficients /,,, 1,,, 12. with 6 are shown 
in Figs. 4,5. For the range of parameters of Table I, /,,, m,, are very small and 
» will be ignored. 


j It will be seen that /,, decreases rapidly as 8—>0 or 1-0, but that it never 
| becomes negative. When h=0-10, w’=0-10, 5=0 or 1:0, /,, is only about 4 per 
cent of the value predicted by quasi-steady theory, but when h=1°0, it is about 
28 per cent of the quasi-steady value. 


Thus, even at a moderate incidence the amplitude of a forced oscillation at 
resonance will be about four times that predicted by quasi-steady theory. 


When 5=~ 0-5, /,, has about the same value as in quasi-steady theory, the 
variations with incidence and reduced frequency being quite small. In fact, only 
when 6 is very small or very close to 1:0 does the influence of w’ and h become 
important. 


As the incidence tends to zero, P and Q —> 0 if 6=0 or 1-0, but /,, does not 
| vanish completely; when o’=0-10, /,,=0-0025 approximately, i.e. the damping at 
zero incidence is almost entirely due to structural hysteresis. On this basis, it would 
| be expected that the amplitude would be largest at zero incidence, but insufficient is 
_ known about the variation of the exciting force with incidence, and it is unlikely that 
| the present analysis is sufficiently exact at small blade angles. No account has been 
| taken of the fact that the vortices must come into contact with the blades—possibly 
' exciting an oscillation—and, if the amplitude of vibration was of the same order of 
' Magnitude as the downward displacement of the wake, non-linear effects would be 
| important. Apart from minor numerical differences, the variation of 1,, with 8, ©’ 
and h is identical with that of /,, and needs no further comment. 


The behaviour of the pitching derivative coefficients m,, and /,, is surprising; 
» m.,, remains virtually constant at 1-0 for all values of 5, ’ and h, thus extending 

one of the results of oscillatory aerofoil theory to rotor blades. But /,. is very much 
| a function of 8, since it changes sign between 5=0-50 and 5=1-0 and, moreover, the 
| peak positive and negative values, although not quite equal, are very large (about 
+20). The peak values decrease with increase of h, and/or of w’, but the sign 
; change is independent of either of these parameters. It is interesting to note that 
| whereas /,, and m,, are both determined largely by the equivalent rotation produced 
| by the wake (see Section 4.2) /,, and /,, are only very slightly influenced by it. 


Thus it has been established that the wake does have a great influence on the 
| damping of flexural vibrations, especially at small angles of incidence, but further 
' confirmation of the theory is necessary, since it is not yet possible to give a complete 
explanation of Meyer’s results and no account has been taken of the shift and 
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Ficure 4. Derivative coefficients—wake theory. Variation of / 1,, with 6, w’ and h. 


12’ 
diffusion of vorticity. Fortunately the vibrations found by Meyer and by the C.A.L. 
were fundamentally different, since the former were forced but the latter were self- 
excited oscillations. 


Hence it should be possible to provide independent confirmation—or 
otherwise—of the present theory, by calculating the flutter speeds of a typical rotor 
blade and comparing the results with experiment. This is the subject of the rest of 
this paper. 

The numerical values of the critical speeds cannot be expected to be accurate, 
since the derivatives are not accurately known, but if the variations shown in 
Figs. 4 and 5 are substantially correct, then the most important experimental 
features should be predicted. 


5. The Flutter of Rotor Blades 
5.1. DERIVATION OF THE FLUTTER EQUATIONS 


The analysis is confined to the binary flutter of a constant chord, symmetrical 
untwisted blade, the degrees of freedom being rigid blade flapping 8 about a hinge 
on the axis of rotation, with centrifugal force providing the restoring moment, and 
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FicurE 5. Derivative coefficients—wake theory. Variation of /,,, with 6, w’ and h. 


tigid blade pitching x about a spanwise axis, coincident with the quarter-chord line 
(see Fig. 6). The elastic stiffness in pitch is supplied partly by the centrifugal force, 
but mainly by a spring at the blade root. This simple model has been chosen for 
two reasons : — 


(i) The experiments of Daughaday and Kline“: © were carried out on a rotor 
of this type. 


(ii) It is a close representation of an actual rotor blade, and if more freedoms 
(e.g., bending and twisting) are taken into account, then the equations of 
motion are very difficult to obtain, and each of the flutter coefficients 
contains a number of small terms, which are functions of the pitch and 
coning angles, and so on. 
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Ficure 6. Blade configuration. 
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The original purpose of the C.A.L. flutter tests was to investigate the influence 
of an elastic coupling between flapping and pitch on the critical speeds. This was 
achieved by attaching a second spring to the blade, ahead of the quarter-chord axis, 
and some distance out along the span,i.e. the stiffness in flapping was partly supplied 
by a spring. If this coupling is included in the analysis, two more variables are 
introduced—the ratio of the distances of the points of attachment from the flapping 
hinge and the quarter-chord axis, and the spring stiffness—and the calculations 
become very laborious. Hence it was decided not to take this coupling into account, 
but to confine attention to the variation of critical speeds with inertia coupling. 
Sufficient tests were done under these conditions to provide a comparison between 
theory and experiment. 


Since the derivation of the flutter equations is standard practice, the details 
are omitted, but some comment on the assumptions is necessary. It is assumed 
that the only significant aerodynamic forces are the lift and the pitching moment 
about the quarter-chord axis. At large angles of incidence this is not quite correct, 
since then the drag forces have an appreciable moment about the flapping hinge but, 
since neither the lift nor the moment are known to any great accuracy, the drag may 
be safely neglected. The aerodynamic inertias and direct aerodynamic stiffnesses 
are all assumed to be negligible, so that the only aerodynamic pitching moment is 
produced by the damping. Since the blades are of high aspect ratio, the pitching 
inertia coupling in the flapping equation is very small and can be neglected. 


Then, putting Y= r, the flutter equations reduce to 


B+nB+B +(4nl../3r) +nljx=0 53) 


/2k?) (8 + X + (n/12k?)*, + =0 
where 410 . . (55) 


and the dots denote differentiation with respect to y. 
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The parameter n, which in simple flapping motion is the ratio of the aero- 
dynamic forces to the inertia forces, is known as the inertia number. 


5.2. SOLUTION OF THE FLUTTER EQUATIONS 


Putting X=x, e” 
(57) 


B = B, ev 
then, for a given o, the unknowns in the characteristic equations of (53) and (54) are 


wo,’ and v. In calculations, however, it is convenient to assume v and to solve for 
and o. 


Then the characteristic equation reduces to 
=mM,{(1 v’)? + n?v7l,,7} / {8 (1 — Gr’ & 


For the usual range of values of the parameters it is found that 1:0 <<» <2°5. 


It will be recalled from Section 2.2 that in the C.A.L. experiments it was found 
that either an elastic or an inertia coupling was necessary to maintain flutter. In 
this case, the possibility of there being an elastic coupling has been ignored, but it 
can be seen immediately that if «=O no flutter is possible. On full-scale rotor 
blades o is usually made as small as is practicable, i.e. the inertia axis is designed 
to lie on the quarter-chord line, but it is of interest to know what is the tolerance 
on o in a particular case. Therefore the critical values of ,;* will be calculated 
for a range of values of c. 


5.3. QUASI-STEADY THEORY 


To estimate the influence of the wake, it will first be assumed that the deriva- 
tives are independent of ’, 6 and h (quasi-steady theory). This is justifiable only 
when w’ is small, but in the present case it is not likely to exceed 0-10 and the 
derivatives can be obtained from two-dimensional theory by writing 


C (iw’)= 1-0, 
where C (iw’) is Theodorsen’s Function. 


Hence it is assumed that 


Mo, and 1,,=0°75. 


The ‘flutter boundaries for the two cases, (a) n=1-414, r’=10-0, (b) n=0-50, 
r=10-0, have been calculated, assuming k?=0-025+0?, and are shown, (marked 
“quasi-steady”) in Figs. 7, 8, 9 and 10. 
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It will be seen from equation (55) that 7? is directly proportional to the spring 
stiffness and inversely proportional to ?. In fact, it is the ratio (natural pitching 
frequency /rotational speed)*, and when presenting results it is convenient to use ©,’ 
since then neither k;, Jy nor 2 need be specified. ©, will be referred to as the 
flutter speed coefficient, an increase in w,” corresponding to a decrease in the critical 
rotational speed. The constant part (1:0) arises from the centrifugal pitching 
moment, which is proportional to 12. 


These boundaries were determined by assuming values for v, and calculating 
@,” and o& by direct substitution in equations (58) and (59). (In the case of “wake” 
derivatives, this procedure is essential, since then /,,, and so on, are functions of ¥), 
The first value of n is typical of full-scale blades. The second is more typical of 
models, and the chosen values of r’ and k*? are midway between model and 
full-scale. 


It will be seen that there is a range of values of o in which flutter is not 
possible, whatever the rotational speed, but the width of this range decreases with 
n. For any given value of o there are two possible critical speeds, the lower one 
corresponding to the case in which pitching is the predominant mode. At the higher 
critical speed, the motion is predominantly flapping, and the frequency ratio lies 
close to that of pure flapping (v~ 1-10). The higher critical frequency shows no 
tendency to remain constant, but increases rapidly with v from about v=1-40, when 
flutter first begins. 


5.4. WakKE THEORY 


The influence of the wake can be estimated if the derivatives are calculated by 
the method of Section 4.2, but the theory is at best qualitative and, as the main 
purpose at present is to establish the characteristics of rotor blade flutter, /,., /.;, 
and so on, need not be known to any great degree of accuracy. Furthermore, it 
has been assumed that the derivative coefficients are constants for the whole blade: 
hence a detailed calculation of the effects of ’ and h is not necessary, provided 
that the influence of 4 is correctly reproduced. Therefore it will be assumed that 
both h and w’ are constants, and it will be seen from Figs. 4 and 5 that the derivative 
coefficients have reasonable mean values when »’=0-10, h=2:0. 


The flutter boundaries for m=0-50, n= /2, r=10-0 are shown in Figs. 7, 8, 
9 and 10; the boundaries predicted by quasi-steady theory are included for 
comparison. 


The most revealing diagrams are those showing the variation of v with ¢ 
(Figs. 7 and 9). Quasi-steady theory predicts a continuous variation, but the “wake- 
theory” shows that there is a range of v in which flutter does not occur. Flutter 
cannot occur for n=0-50 in the range 16<»< 2:0, and for n=1-414 when 
1:8<v<2-0. Flutter may be possible in these frequency ranges if o is negative, 
but calculation shows that the centre of gravity must be well ahead of the quarter- 
chord line in the present case. Thus, the most important experimental result that 
two different flutter modes are possible at the same c.g. offset is amply confirmed. 
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For convenience, flutter when v > 2:0 will be referred to as the high frequency 
or h.f. mode and flutter when 1-0 < v < 1:8 will be referred to as the low frequency 
or lf. mode. The minimum value of o necessary to cause flutter will be 
denoted by on. 


It will be seen that o,, is different for the two modes, the difference depending 
on. (All four values are different from those predicted by quasi-steady theory, 
but the difference between the quasi-steady and 1.f. results is quite small. This is 
because /,, 1-0, and /,, does not have much influence on the results at these 
frequency ratios, see equation (58). When n=0-50, o,, (hf.) is about 0-06, and 
m (1.f.) is about 0:04, but when n=1-414 the two values of o,, are almost equal. 
Thus, if the range of c.g. travel in an experiment is limited, the high frequency mode 
may not be excited if n is small, but it may be found at larger values of n. This 
point is confirmed by the C.A.L. experimental results. On the rigid (i.e. stiff against 
bending) blade which had a very small value of n, the high frequency mode was not 
excited over the available c.g. range, but on the flexible blades this mode was 
produced even at 4° incidence. These blades had a value of n ~ 0°75, against 0-25 
for the rigid blade. Thus, differences in the flutter characteristics of the two sets of 
blades cannot, as was assumed by the authors*:®, be attributed directly to 
flexibility. There can, however, be no doubt that flexibility is important at 
frequency ratios greater than 2:0, because this is close to a natural bending 
frequency. 


The variation of the critical speeds with o is quite complex, since for a given o 
there may be anything up to four speeds. For n=0-S0, the critical speed boundaries 
of the two modes intersect at a certain value of 7, so that, when o exceeds this, two 
critical speeds, the lower of the two for the h.f. mode and the higher for the 1.f. mode 
speeds, will be obtained. If flutter is excited at this value of o, three critical speeds 
will be obtained, one of which contains both modes, and therefore frequencies, 
together. For o less than this, there may be four, three, two, one or no critical 
speeds. When n=1-414, the range of o in which four speeds are possible is much 
smaller, and beyond the intersection of the two boundaries there is not much 
difference between the lower critical speeds of the two modes. Also the range of 
speeds in which the high frequency mode can occur is reduced. 


6. Concluding Remarks 


6.1. GENERAL 

There is sufficient agreement between the foregoing results and experiment, and 
between the calculations on damping and Meyer’s experiments, for it to be claimed 
that the wake is primarily responsible for the observed phenomena, but the analysis 
is necessarily limited and there are still many points to be explained. It is probably 
beyond the scope of any reasonable theory to take into account all the various 
factors (incidence, Reynolds number, diffusion of vorticity, number of blades, and 
so on); only controlled experiment can do that. 


The predicted variations of m,, (constant), /,. (sign change) and /,, could easily 
be detected in straightforward experiments, but the tests would have to be done in a 
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compressed air tunnel, to keep the Reynolds number constant, and it would be 
necessary to provide for an independent variation of stiffness to keep w’ constant. 


Implicit in the theory of Section 3.2 is the assumption that the conditions which 
exist when N<v=<=N+1, are the same for all values of N, but it will be 
remembered that, whereas Meyer found vibrations for values of v as large as 8:0, 
Daughaday and Kline obtained “wake flutter” only for v=2-0 approximately. How 
much of this was due to the different recording techniques is not known, but the 
foregoing assumption is not completely justified. The larger values of v are obtained 
only at low rotational speeds, when the total energy available is small, the Reynolds 
number is small, and the “critical” speeds are very close together. Thus it may be 
difficult to excite the oscillations at large v unless there is a strong exciting force. 


6.2. AIRSCREW FLUTTER 


The emphasis in the preceding sections has been on the rotor blade, but there 
is no reason why the wake should not influence the flutter and vibration of airscrews 
to some extent. The theory of Section 3.1 would not be valid, since the aspect ratio 
of airscrews is not as high and, except on the spinning tower, the pitch of the helix 
is quite large. So far no cases of “wake” flutter on airscrews have been reported 
as such, but the following quotations from Lynam’s® early paper describe a 
vibration so similar to those of Ref. 3 that it can only be concluded that flutter of 
this type is possible. 


“A certain airscrew .. . . flutters in a rather unusual manner.” 


“It has definite critical speeds at which flutter occurs and between which 
speeds the flutter is not so noticeable, suggesting that in this case the 
fundamental is below the running speed.” 


“A special airscrew, having no pitch and symmetrical sections of the same 
thickness, at the same proportions of the chord, and also having the same 
variation of chord with radius, was constructed and spun.” 


“Although, owing to the pitch being zero, and the sections symmetrical, 
the lift force was eliminated and there was no thrust, the blades fluttered at 
practically the same speeds as the airscrew.” 


Lynam states that the method of eliminating the flutter was to increase the 
thickness of the sections over the outer third of the blade. This had the effect of 
greatly increasing the natural bending (as well as torsional) frequencies, especially 
in rotation, and thus removed the critical speeds outside the working range. 


6.3. THE ROTOR IN FORWARD FLIGHT 


Finally, it is necessary to consider the influence of the wake on the blade motion 
of a helicopter in forward flight. Once again the theory of Section 3.1 is inapplic- 
able, because the wake now moves rearwards as well as downwards, and the 
relative velocity is periodic, but it seems likely that the wake still has some 
influence. 
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FLUTTER OF ROTOR BLADES 


The experiments showed that the vibratory stresses have a peak amplitude at 
low forward speeds, i.e. in the “transition” region, when the induced flow is not 
uniform around the disc. Thus an oscillation will be forced at low », but it is 
unlikely that the effects of the wake will vanish abruptly as the forward speed 
increases from zero. The conditions for minimum damping will no longer be 
realised exactly, but there will be some decrease in /,, at integral values of v, and 
this, together with the increased forcing function, will cause an increase in the blade 
stresses. As the forward speed increases, the forcing function will at first decrease 
—after leaving the “transition” region—but will later increase rapidly. Presumably 
the aerodynamic damping will increase with forward speed, but the remaining 
blades and their wakes may play an important part in determining the aerodynamic 
forces, and it is not possible to draw any further conclusions without additional 
tneoretical and/or experimental evidence. 

All previous published attempts at the calculation of the forces on a blade in 
forward flight have assumed that the blades were in steady motion. It is now 
apparent that this assumption must be dropped, if reliable estimates of the stresses 
are to be obtained, but the labour involved would be enormous and the most 
satisfactory method would be to measure the derivative coefficients in a 


wind tunnel. 
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Appendix 


+00 
EVALUATION OF THE INTEGRAL + n7h?). 


Integrating by parts gives 


+00 +00 
[ / (2 + | + +n?h?)] dn] io’ 


+00 
=0+ | + 2n?) (4? + n7h?)?) 
+00 +90 
=< | (n? + n2h?) +2h?n? /(n? + n*h?)?, 
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These integrals are standard forms and may be obtained from tables of definite integrals, 


+00 

+00 

j edn +n2h?)? + 20’nh} /(2nh)? 

+00 

(n? + n2h?) = e 
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The Lift on a Thin Aerofoil with a 
Jet-Augmented Flap 


D. A. SPENCE 


(Royal Aircraft Establishment) 


SUMMARY: The method which was developed in an earlier paper") for calcu- 
lating the lift on a thin jet-flapped aerofoil has been applied to the case in 
which the jet enters the main stream tangentially from a deflected flap. 
The lift derivative with respect to the (small) angle of deflection 2 is of the form 


0C,/0B=2(x +sin x +2zD,), 


where y=2 sin~',/E, E being the ratio of flap length to aerofoil chord, and 
D,, a function both of E and of the jet momentum coefficient C,. In the special 
case C;=0, D, vanishes and the classical results for attached flow on flapped 
aerofoils without super-circulation are obtained. More generally, D, has been 
evaluated over the range 0<(C,<5 for six ratios E between 0-02 and 0-50, 
using an electronic computer. The lift derivatives so found lie between those 
in the limiting cases of “deflection” (E=0) and “incidence” (E=1) which were 
calculated in Reference 1. They show that the effectiveness of the jet-flap 
scheme might be greatly increased by blowing over flaps of the order of 20 or 
30 per cent of the chord, rather than over the small “shrouds” which have 
hitherto been proposed. 


l. Introduction 


Improvements in the lifting effectiveness of the jet-flap scheme could arise if 
the air were ejected over a trailing-edge flap of moderate size (flap blowing), instead 
of from the trailing-edge (pure jet flap), or over a small flap merely intended to 
turn the jet (shrouded jet flap). The author’s treatment of pure jet flaps by thin 
aerofoil theory” is here extended to flap blowing, and comparative results are 
derived for various ratios of flap to aerofoil chord. The practical significance is 
discussed in Section 5 mainly in relation to jet-flap schemes, where the momentum 
ejected over the trailing-edge flap exceeds that required to prevent flow separation 
by a large amount; but brief reference is also made to more conventional boundary- 


layer control applications, where the excess momentum leaving the trailing edge 
is small. 


In the earlier theoretical treatment™’, an inviscid flow solution was obtained 
for the case of a thin deflected jet of momentum flux coefficient C,, emerging at a 
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small angle 8 to the chord, from the trailing edge of a thin aerofoil whose surface 
slope is prescribed by 
dy 


1. 


For the case of flap blowing, we now write 
e(x)=0, O< x<(1-E) 
e(x)=B8, (1-E)V<x<l 


and allow the jet to emerge tangentially at the trailing edge. 


For C,=0, the Glauert™ theory for hinged flaps is reproduced, but the results 
are obtained in a more compact form, as a result of using the Carlemann™ inversion 
formula for Cauchy integrals. As this treatment appears to have some intrinsic 
interest, and also provides a simple illustration of the general method, it is presen- 
ted first in Section 2. The jet case is then presented in Section 3. Results have 
been computed on the Deuce electronic machine for a range of C,; values and a 
number of flap/chord ratios. These are summarised in Section 4 and discussed 
in Section 5. 


NOTATION 
The chord length is taken as unity throughout. 


x,y co-ordinates of aerofoil and jet sheet in directions parallel and 
perpendicular to the undisturbed stream 


£, dummy variables over the same range as x 
py.U, density and velocity of undisturbed stream 
w(x) downwash velocity on y=0 


f(x) non-dimensional vorticity distribution on <1 representing 
wing (i.e. velocities are U,{1+4f(x)} on upper and lower surfaces) 


jet momentum coefficient, =J/(4p,U,7) where J=flux of 
momentum 


A= 4/C; 


g (x)= —(C,/2) (x), non-dimensional vorticity distribution on x > 1 
representing jet sheet 


¢(x) slope of aerofoil surface 
a incidence of forward part of surface 


8 angular deflection of flap from forward part of surface 


E ratio of lengths of flap and whole aerofoil 
6 cos-'(2x-—1) for 0<x=< 1 

x value of @ at hinge point, i.e. cos~' (1—2E) 
/{1+(1 - 
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¢=2cos-'(x-#) for x>1 
k(x) defined by equation (20) 
K (¢)=[1/(28)] g [1/cos? (¢/2)] 
number of pivotal points 
dm=mt/N, m=0,1,...(N—1) 
.Dy-, coefficients in Fourier series for K’ (¢)/sin (¢/2) 
Anns Omns fm coefficients in equation (26) 
C, lift coefficient 
Cy nose-down pitching moment coefficient about leading edge 


Cg moment coefficient of flap about hinge line 


2, Attached Flow on a Thin Aerofoil with Hinged Flap (No Jet) 


Consider an aerofoil on the real (x) axis from 0 to 1, with a hinge at the point 
x=1—E, where it is bent through a small angle ~. 


X=0 (I-E) 


FiGcure 1. 


Suppose the aerofoil is represented by a vortex distribution of strength U, f (x) 
per unit length, where U, is the velocity at infinity. 


The downwash on the axis is then given by 


On the segment (0, 1) this is known to be 
w (x)=0, 0<x<(1-E) 
w (x)=U,8, (l-E)<x< 


Inverting equation (1) by means of Carlemann’s identity* gives 


*The general solution of equation (1) would include a term C/{x(1—x)!} on the right-hand 
side of (4), C being an arbitrary constant. But, in accordance with the Joukowski hypothesis, 
we take C=0. 
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x 


Substituting from (2) and (3) 


4 


Writing x=4(1+cos 4). 
1—-E=4(1+cos x), 
the Cauchy integral may readily be evaluated, and the loading is given by 


f(x) _ 9 ,.1, sin (6+ x)/2 


The total lift coefficient is given by 


1 


C,=2 { f (x) dx. 


0 
It is found by integrating equation (8) that 


(This equation may also be obtained directly by integrating equation (5). Note 
that with x=7, i.e. E=1, it gives the usual incidence result.) 


It may readily be verified that Glauert’s formula for lift is identical with (10), 
the necessary changes in notation having been made. In addition, it may be shown 
that the Fourier series for the loading given by Glauert has the sum given by 
equation (8). The various moment coefficients, on the aerofoil and flap, may also 
be obtained by integration of equation (8), giving the following results. 


The lift coefficient on the flap (referred to the aerofoil chord) is 
Ci tap = 28x?/7.. 
The nose-down pitching moment about the leading edge is 
Cu=B8 (4x.+ sin x sin 2y). 


The contribution to this pitching moment from the flap is 


Cy flap — (x + sin x). (13) 
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JET-AUGMENTED FLAP 


The hinge moment on the flap is 


Cu Cy flap — = E) Ci flap 


(14) 


(Here the hinge moment is referred to the aerofoil chord; to make it non- 
dimensional in the usual way one would divide by E’, =sin* (4x). 


Similarly, in linear theory the effect of incidence is directly additive; thus, with 
an incidence z the lift coefficient becomes 


3. Flap with Tangential Jet 
3.1. INTEGRAL EQUATION FOR SLOPE OF JET 


The case when a jet is blown tangentially over the flap from the hinge line may 
be treated precisely as if the jet emerged tangentially at the trailing edge, provided 
that an effectively attached flow is obtained. The effect of the jet on flow over the 
flap itself is merely to prevent separation, and the pressure distribution is not 
modified in any other way. (The pressure rise across the jet is proportional to its 
curvature, and thus is zero for the part in contact with the flap.) 


The analysis of Ref. 1 is thus applicable to the system. In that paper the 
deflection of the jet below the axis is taken as y(x), (1 << x <0), and the local 
gradient of the aerofoil surface is :(x), (Q<x<1). The jet momentum-flux 
coefficient 


C,=J/(4p,U,"). ‘ ‘ (16) 
is defined. where J is the momentum flux per second entering the jet, and p,. U, 
tefer to the undisturbed stream. 
With the transformation 


A=4/C, 
(17) 


2 
y’ (x) 


the following singular integral equation is obtained for the slope function g(x) 
(equation (69) of Ref. 1) valid in 1 < x << 00:— 
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If < (x) is defined by the step function 


(x)= | 


8, cost X<x<i| 


the right hand side of (18) is found to have the value 28k (x), say, where 


so that equation (18) is 


where k(x) is known and g(x) unknown. 


3.2. SOLUTION OF EQuaTION (21) 


This equation is of precisely the same form as that solved in Section 5 of 
Ref. 1; it may be transformed by writing 


x=1/cos* 


(x)= 28 K (9), 


where K’ (#)=sin 5 by D,, cos n¢. 
n=0 


With this transformation, (21) becomes 


+ 3 D,, cos* sin nd +AK (9) = —k [1/cos? (@/2)]. (25) 
n=1 


D, cos? sin 
By satisfying this equation at N pivotal points 


mr 
m=0, 1,2... (N-1), 


N linear equations are obtained for D, ... Dy_,, of the form 
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1 
+ ABmn) Ds = fu m=0, 1, ... (N-1), . 


where Ang =SIN  Ann=(1+Cos sin nd, for n>O 


(cos no», +2n tan sin ) 


as in Ref. 1, and 
fu= —2 sec [1 /cos* (%m /2)] 


= Bie ® Pm 


The set of equations (26) with N=9 has been solved for a number of values of 
C; between 0:05 and 5 for flaps of flap/chord ratio E (i.e. sin® (4x)) equal to 0-02, 
0:05, 0:10, 0:20, 0:30 and 0-50 on an electronic computer. The set of coefficients 
{D,,} converges rapidly in each case; e.g. for C;=2, E=0°3, D,, D, ... Ds =0°3513, 
0:3013, 0-:0008, 0:0270, 0:0011, 0:0066. — 0-0006. 0-0026, —0-0005. 


3.3. EXPRESSIONS FOR LIFT COEFFICIENT AND CHORDWISE LOADING 


The chordwise distribution of vorticity is taken as U, f(x), (0< x <1), as in 
Section 2, and from equation (62) of Ref. | it is known that 


The lift coefficient is given by 
1 


f()dx+BC,, 


the term BC, representing the direct vertical momentum component of the jet. The 
contribution to the integral of the term involving the step-function <() is already 
known; as in Section 2, equation (10), it is 28(x+sin x). The contribution from 
the Fourier series is evaluated as in Section 6.1 of Ref. 1 and it is found that 
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TABLE I 
COMPUTED VALUES OF 


0-05 | O2 0:3 | 10 


| 1:774 | 2-487 3-455 | 4151 | 5: 6283 
1-845 2:547 3-505 | 6°392 
| 1: 2-081 2:751 | 3-681 4:360 | 6°460 
2:318 | 2-962 3-868 4:535 | 6614 
1641 | 2: 2-708 | 4-190 4839 | 6-885 
2:391 3-340 3-908 | 4740 | 5-365 7:365 
2:707 | 3 3-618 4173 4-989 5-607 7-588 
4-026 4:827 5-337 6:102 7:561 8-605 
6836  =7°301 8-012 8-566 10°405 
11-231 11-825 | 12:236 12°878 13-385 | 14-227 15-109 


Note: The values for E=0 and E=1 are those of 0C,/07 and OC, 182 listed in Table 3 of 
Reference 1. Very close interpolations to them are provided by the formulae 


AC,,/d7 [1 +0-151 +0139 C,]} 
C,3 +0:219C, 


Thus the lift coefficient is given altogether by the simple formula 


(x +sin x+2zD,) B, 


where y=2 sin-'/E, E being the flap/chord ratio, and D, depends on E and C,. 


In the same way the chordwise loading is made up of a closed term coming 
from the solution for a Glauert type flap, equation (8), together with a series 
identical to that obtained in Section 6.2 of Ref. 1; in fact 


where 


3.4. EFFECT OF OVERALL INCIDENCE 


Since this is a linearised theory, the effect of an overall incidence z is directly 
additive; thus the expression (32) becomes 


where B, is the coefficient evaluated as a function of C, in Ref. 1. It may be 
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FiGurE 2. QC,/08 for thin wings with blown flaps. 


verified that the same result is obtained by taking the step function < (x) as 
(x)=2, 
=a+B, cos? x 


Similarly the more general expression for chordwise loading is derived from 
(33) by adding 22 {(1—x)/x}! to f (x) and replacing 28D, by 2 (2B, + 8D,). 


4. Results of Calculations 


The computed curves of 0C,,/08 which, by (32), is given by 


oC, 

3B 
for E=0-02, 0-05, 0-10, 0:20, 0:30 and 0-50 are plotted against C, in Fig. 2, and 
their numerical values are given in Table I. The figure shows also the curves of 
0C,,/8a and 8C,/ 27 obtained in Ref. 1. All the new curves lie between these two, 
and it is clear that 0C,,/d7 may be regarded as the limit of 0C,/08 as E—> 0, and 
0C,,/02 as the limit when E —> 1. 
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FiGurE 3. Nose-down pitching moment about leading edge on aerofoil with blown flap. 


2 3 4 5 


Ficure 4. Lift and hinge moment coefficients on 10 per cent and 30 per cent blown flaps. 
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Ficure 5. Blowing over a 30 per cent flap. Vorticity distribution. 


Ficure 6. Blowing over a 30 per cent flap. Non-dimensional jet shape. 


Figures 3 and 4 show the pitching moment and hinge moment curves for 
E=0-10 and E=0-30. These have been obtained by integrating the contribution 
of the Fourier series to the loading expression (33) numerically and using the 
results in Section 2 for the contributions not arising from the jet. (The sum of the 
Series is non-singular, unlike the “Glauert” term, and the numerical integration is 
straightforward; for the pitching moment the integration can in any case be 
performed analytically, as in Ref. 1). 


Figures 5 and 6 show the distribution of vorticity on the aerofoil and jet, and the 
jet shapes respectively, for a number of C, values and E=0-30, It is seen that the 
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vorticity distribution is singular at the hinge point. The singularity* comes only 
from the “Glauert” term 


ios + (=z) | 


(3) 


5. Discussion 


Comment should be made about the applicability of the analysis of Ref. | to | 
the present scheme. The C, entering the analysis is the effective jet momentum 
flux at the trailing edge; this will obviously be rather less than the momentum 
emerging at the slot at the hinge point, since some momentum is transferred to the 
boundary layer to prevent it from separating in the adverse pressure gradient on 
the flap. This momentum loss is determined by the internal mechanics of the jet 
flow over the flap and is beyond the scope of the present treatment, which is 
confined to the external inviscid flow, the physical properties of the jet having 
entered the calculation solely by means of the value xU,C,/2 of the vortex strength 
on 1<x< © used in forming equation (18) («=curvature of sheet ~ d*’y/dx’). 


The present assumptions are entirely consistent with those underlying the 
treatment of the pure jet flap problem in Ref. 1. Assuming an infinitely thin jet 
in an inviscid flow, it may be shown that the variation of the momentum flux 
along the jet tends to zero, however the static pressure may vary. In this limiting 
case, therefore, there is no loss of jet momentum between the slot and the trailing 
edge, and the present results are strictly true. 


In practice a certain “threshold” momentum must be necessary to establish the 
prescribed flow in either case, and this may well be greater in the case of the blown 
flap than for the pure jet flap itself. However, it would appear, in general, that it 
is the effective momentum flux at the trailing edge that determines the “super- 
circulation”, and hence that the present results are relevant to a real flow provided 
that 


(i) the prescribed “attached” flow is attained in the real flow, and 


(ii) the relevant jet momentum flux in the real flow is taken to be that at the 
trailing edge. 

If the “true” C,’s are defined at the slot, they must be greater than those shown 

in Fig. 1; how much greater would depend upon the development of the jet over the 


*In fact on a plate bent through an angle 8 the singularity should be like that of 
(x—x,)~A/ 


at x=x,, but the assumption of a step-wise discontinuity in downwash produces instead the 
logarithmic singularity of (38). 
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flap in any particular problem. For zero “true” C,, the prescribed flow is simply 
not attained, because of flow separation, and the real values of 0C;,/08 would be 
substantially lower than the intercepts shown for C,=0. 


Thus certain unknown factors limit the conclusions that may be drawn from 
the results plotted in Fig. 1, but it appears there that, by the use of a blown flap, 
either 


(i) with a given C, and deflection, higher C,’s are obtained with larger flaps; 
e.g. for C,;=1, 50 per cent more lift is obtained with a 20 per cent flap 
than with the pure jet flap; 


for a given C, and deflection, the larger the flap, the smaller is the 
required C,; e.g. for 8B=4 radian=28-6°, Cy=3 could be obtained at 
C,=2 with a pure jet flap, or at C;=1 with a 20 per cent chord blown flap. 


These are substantial gains which may well be attained, in part at least, in a 
real flow, provided that an efficient slot-flap combination can be designed. Thus 
it appears that the effectiveness of the jet flap scheme might be greatly increased 
by the use of flaps of finite size (substantially greater than the “shrouds” which 
have been discussed). At the opposite extreme of low C,’s, it seems that useful 
gains might be obtained from “super-circulation” in the conventional boundary- 
layer-control applications. 


It may be noted that the 2 per cent flap is sufficiently different from the pure 
jet flap to indicate that trailing edge conditions can make a good deal of difference. 
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Some Identities on Structural Flexibility 
after Buckling 


E. H. MANSFIELD 
(Royal Aircraft Establishment) 


SuMmMaRY: This paper gives some identities concerning the flexibility of certain 

elastic structures, including flat plates, at the onset of buckling and subsequently. 

The identities have a useful application in determining the flexibility and the 

rate of change of flexibility with load for structures in which the mode in the 
buckled state changes continuously with load. 


1. Introduction 


A flat plate, on being loaded in its own plane, deforms initially in its own plane 
and possesses a certain flexibility. As the loading is increased, the plate may 
buckle out of its plane, and there will be a sudden increase in the flexibility at the 
onset of buckling. As the loading is further increased the out-of-plane deflections 
will increase, but if the out-of-plane mode does not alter then the flexibility will not 
alter, and the load-deflection and load-flexibility curves for the structure will be 
as in Fig. 1. 


If the out-of-plane mode alters continuously with the loading, as the loading 
is increased beyond the initial buckling stage, then the flexibility will alter, and 
the load-deflection and load-flexibility curves for the structure will be as in Fig. 2. 


FLEXIBILITY 


DEFLECTION 


IF NO BUCKLING 


LOAD LOAD 


Ficure 1. Load-deflection and load flexibility curves for structure with constant mode after 
buckling. 
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DEFLECTION FLEXIBILITY 


LOAD LOAD 


FiGuRE 2. Load-deflection and load-flexibility curves for structure with varying mode after 
buckling. 


The calculation of the flexibility for the structure considered in Fig. 2 is 
considerably more difficult than for that considered in Fig. 1, and it is the object 
in this paper to derive certain identities which enable the behaviour of such a 
“varying-mode” structure to be expressed in terms of the behaviour of a “constant- 
mode” structure. The identities were useful in an earlier paper by the author 
in connection with the post-buckling behaviour of reinforced plane sheet under 
shear. 


NOTATION 
A,r’ applied load 
deflection in the direction of A 
flexibility, =dA/daA 
work done by the applied load 
flexibility before buckling 
actual buckling load 


buckling load for the structure if it were constrained to buckle in the 
mode which the structure actually would have at a load A(A > JA,) 


flexibility of the structure constrained to buckle in the mode which the 
structure actually would have at a load A 


deflection not in the direction of A 
“ cross-flexibility,” —dA’/dA 


Suffix “cc” refers to conditions at the onset of actual buckling. 


2. Derivation of the Identities 


2.1. GENERAL CONSIDERATIONS 


Consider a structure subjected to a load 4’, greater than A., and suppose that 
the out-of-plane mode has altered from M., the initial buckling mode, to M. The 
state of the structure at this load is independent of the loading history and would 
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“ 

Xn x Am 

FiGurE 3. Load-deflection and load-flexibility curves for structure in constrained 


and unconstrained state. 


be reached if there had been present in the structure constraints which prevented all 
out-of-plane modes except M. Such a constrained structure would buckle at a 
load 4,,, say, and would have a constant flexibility 9,,, say, after buckling. For 
A <4, there would be no difference between the unconstrained and the constrained 
structure. The unconstrained and constrained loading paths are indicated by the 
full and broken lines in Fig. 3. 


The fact that the state of the structure is independent of the loading history 
can be expressed by the equations 


and . . . . (2) 
since the internal energy stored may be equated to the work done in each case. 


Furthermore, 4 and T can be expressed in terms of \’ and ¢:— 


A 


A= dr, 


T= dX. 


Now, in the constrained state, 


$= for A<Am (4) 
for A>A,,, 
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so that equations (1) and (2) may be combined with (3) and (4) to give 
| dA= PyAm + Pm Am) 
0 


0 


Equations (5) and (6) are valid for any value of ’, and may therefore be 
differentiated with respect to A’ to give expressions for »; furthermore, now that 
\’ is the independent variable it may be written for convenience as A. 


The main identities are to be derived by differentiating and combining 
equations (5), (6) and using the condition that, at the onset of buckling, 


2.2. THE IDENTITIES 


2.2.1. Identities in the Buckled State 


dd, 


=Om + A - An 


do» dx 


Identities at the Onset of Buckling 


Identity (C) means that the flexibility at the onset of buckling is independent of 
the rate of change of the shape of the mode. 


Identity (D) is a restatement of the fact that at the onset of buckling the true 
buckling mode is such as to make the buckling load a minimum’. For the same 
reason (d°A,,/dA*),. cannot be negative. 
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Since . >, and (d°A,,/dA*). is not negative, it follows from identities (E) 
and (F) that (d¢/dA). is not negative. This means that the flexibility immediately 
after the onset of buckling cannot decrease with increase of load. 


2.2.3. Identities for Cross-Flexibilities 


Similar identities may be found if A’ is a deflection in a direction different 
from A, but such identities are now to be deduced only from equation (5), as equa- 
tion (6) was based on equating work done. 


Identities corresponding to (C) and (E) are now 


pe (Om). 


3. Conclusions 


Certain identities have been derived in which the post-buckling flexibility of a 
structure with a post-buckling mode that varies with the load is related to the 
post-buckling flexibility of a structure constrained to deflect in a constant post- 
buckling mode. 


The identities have a useful application in determining the flexibility and the 
rate of change of flexibility with load at the onset of buckling. 
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A Note on the Drag of Wing-Body 
Combinations in Linearised Supersonic Flow 


L. E. FRAENKEL 


(Department of Aeronautics, Imperial College) 


Summary : In this note Ward’s ‘“area-transfer rule’ is modified in order to make 
it more suitable for configurations with fuselages or nacelles which are neither 
slender nor smooth (although they are assumed to have small profile slope). 
The body source strength of Ward’s analysis, which is appropriate only for 
smooth, slender bodies, is replaced by a source strength which is valid more 
generally. However, the drag contributions of the body multi-sources and of 
the interference field are still neglected, so that the theory is tentative rather 
than rigorous. 


Introduction 


Substantial advances towards calculation of the wave drag of whole aircraft 
configurations have recently been made as a result of work by Hayes“ and others 
in the United States, and by Ward in the United Kingdom. The present note deals 
with a slight modification of Ward’s principal result, the so-called “area-transfer 
tule.” 


A wing-body combination free of local transverse forces may be represented 
by a surface distribution of sources and a line distribution of sources and multi- 
sources: other representations, of course, are also possible. Two sets of line 
singularities are required in general: those appropriate to the body alone, which will 
be called “body” sources and multi-sources, and those required to cancel the normal 
velocity induced on the body by the wing, which will be called “interference” 
sources and multi-sources. The strength of the wing sources is immediately known, 
but, in general, the strength of the line singularities, and their contribution to the 
drag, are not calculated so easily. However, if the body is “slender” and “smooth” 
considerable simplifications result: these form the basis of Ward’s analysis. 


A precise explanation of terms is now required. By a slender body we mean 
one whose maximum radius is small compared with both its length and some 
characteristic wing dimension, like the chord or the square root of the plan form 
area. Throughout this note, bodies with non-vanishing nose and base areas are 
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taken to be extended to infinity upstream and downstream by cylinders parallel to 
the free stream: a smooth body is then one such that the slope of the resulting 
profiles is continuous everywhere. Only “net” wings will be considered; that is, no 
hypothetical extension of the wing inside the body will be contemplated. 


The simplifications already mentioned are as follows. If the body is slender, 
then 


(i) the drag contribution of the interference sources and multi-sources is 
negligible, 


(ii) the drag due to interaction of the body multi-sources with the wing 
sources is negligible. 


If in addition the body is smooth, then 


(iii) the body source strength is simply the rate of change of cross-sectional 
area, 


(iv) the drag contribution of the body multi-sources is negligible. 


The drag therefore consists of only three terms: the drag of the wing alone, the 
drag of the body sources alone, and a coupling term between the wing sources and 
the body sources. 


Most aircraft fuselages are neither slender nor smooth; on the other hand an 
“exact” analysis of the drag (that is, one which is exact within the framework of 
linearised theory) leads to formulae which are not suitable for computing the drag 
of general configurations. (Such an exact analysis is made in Ref. 3 for wing-quasi- 
cylinder combinations). A compromise is therefore sought in the present note. The 
body source strength of Ward’s analysis, which is appropriate only to smooth 
slender bodies, is replaced by a source strength which is valid more generally, but 
it is assumed that the statements (i), (ii) and (iv) are still approximately true. The 
new body source strength is calculated by combining the quasi-cylinder and slender- 
body source strengths in the manner already used for the pressure on the body by 
Lighthill (Ref. 4, p. 454), Ward (Ref. 5, p. 216) and Warren and Fraenkel, 


The drag result of this note requires more computation than does Ward’s 
original form, because the line source strength must be computed first, but this may 
not be too difficult if an electronic machine is used. A more serious defect of the 
new drag integral is that it is much more difficult to minimise than Ward’s 
expression, because it involves the line source strength rather than the body 
geometry directly. The question of whether these new complications are worth- 
while can probably be decided only by trial. 


NOTATION 
B=(M? - 1} 
drag force 
F(z) special function, see equations (3) to (5) 
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body source strength 

Heaviside unit function 

transferred area function, see equation (9) 

modified Bessel functions of the first and second kinds 
Heaviside operator with respect to z 

free stream Mach number 

body radius 

position vector 

vector distance from z-axis 

cylindrical co-ordinates with z-axis parallel to free stream 
body cross-sectional area 

wing thickness 

free stream velocity 

free stream density 


perturbation velocity potential 


2. The Strength of a Line of Sources 


A fundamental problem of the flow past bodies of revolution at zero incidence 
is: given the potential 


z-Br 
U j g (z,) 


(r, Z)= - Vi(z-z,)? - Br] dz,, 


co 


to find the source strength g (z,) such that the boundary condition 


op _ 
ap (R @). = UR’ (2) 


is satisfied. 


With a quasi-cylinder, for which R (z) is nearly constant and equal to R,, say, 
the solution of this problem is as follows. We introduce the function* F (z), 
defined by 
=F, (2) H (2)=F (2) (3) 
pK, (p) 


where H is the Heaviside unit function. 


*The derivative of this function has been roughly tabulated by Nielsen‘*. 9, who denotes our 
F(z) by —M,(z). Since the integral in relation (8b) can be written in terms of F’ and S’, 
instead of F and S”, F’ may, in fact, be the more useful function for practical purposes. 
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The behaviour of this function is indicated by 


F (z)~ 160 


| 1 


29 
~1- 5 t E ~ 3 log 2(2-1)| 


and the function could be tabulated from the identity 


co 


F (z)=1- fere-n 1, Q) dx 


Then the quasi-cylinder source strength is 


(Z= F (25 + 1) dS’ (z,). 


The source strength of a smooth slender body is 


Both these results are approximately contained in each of the following, which can 
be shown rigorously to be valid approximations for any body of revolution with 
continuous radius and small profile slope : — 


| F +1) as" (eo. 


—CS 


| F 


g. and g,, which are only two of a large number of possible functions, are 
equivalent because F differs significantly from unity only when its argument is not 
large, that is when z-z, is O(R) or less, and in these regions R (z) =~ R (z,). It 
might appear at first sight that g, should correspond exactly to the Lighthill form 
of the pressure on the body, while g, should correspond to the Ward and Warren- 
Fraenkel form, but this is not the case; further justifiable approximations are 
required before either of the source strengths given yields precisely either of the 
forms for the pressure on the body given by these authors. 
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It is important to notice that, in the region of the cylindrical extension behind 
a body of finite base area, the source strength does not vanish: sources are required 
in this region to make the flow parallel to the free stream, but their strength can be 
shown to decay like z~? for large z. Behind a pointed tail g, vanishes exactly, but g. 
only approximately; g, could, of course, be “cut off” in such a region. 


The source strengths g, and g, are continuous: a discontinuity in S’ (z,) at 
z,=e; contributes a term which behaves initially like 


and after a few radii like AS/. 


3. The Representation of a Wing-Body Combination by Singularities 


In what follows it is assumed that the mean wing surface and the body axis are 
both parallel to the free stream. 


The wing may be represented by sources of strength T’ (R) distributed over its 
mean surface. Alternatively, it may be represented by line sources and multi- 
sources distributed along the body axis: these give a valid representation of the 
flow outside a circular cylinder which is centred on the body axis and contains the 
wing. The line source strength in this case is 


rm) 


| z- | < 8r, 


which is Ward’s transferred-area function and is denoted by A’ (z) in Ref. 2. 


The former representation of the wing is useful for calculating the drag of the 
wing alone; the latter is useful for calculating the drag due to interaction of the wing 
and body fields. 


The body is represented by a line of sources which satisfies the body boundary 
condition averaged around the body contour: the body multi-sources and the inter- 
ference sources and multi-sources are neglected. If R,,(z) is a mean radius, 
averaged in some convenient way with respect to 6, the body source strength is 
taken as 


| F 1) ds’ 


or as the corresponding g (2). 
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4. The Drag 


Upon application of the methods of Refs. 1, 2 or 3 to the present problem the 
drag is found to be 


T” (R,) T” (R.) cosh- + 


(z,) g (22) log | z, | dz,dz. + (z,) (2,) log | z, - z. | dz,dz. | 
(il) 


| T” (R,) T” (R,) cosh- <2 dS dS, + 


+ (z,) +H (2a) +H log | 2, | - 


- {J h’ (z,) h’ (z,) log | - |dz,dz. 


A new feature of these results is that the third integrand of (11) and the second 
integrand of (12) vanish only outside semi-infinite (instead of outside finite) regions. 
But the integrals can be expressed as finite integrals by a simple change of variable. 
If the origin is chosen so that g(z) and /(z) vanish for z< 1, these integrals are 
of the form 


| f’ log | z, | dz,dz,. 


Putting z,=1,~', z,=t.~', we have 


1 1 1 


1 
1 “a d d 
I= =| dt, f (t,~ dt, f (t.~') log |t, - t,|dt,dt. - f(t, | ) dt. - 
0 0 


1 


1 
d d 
"Mog | | dts | )dt, 
0 


0 


since g (z) and A (z) vanish at z=1 and z=00. (d/dt)g(t-') is O(t) for t—> 0 and 
O (1 -1)-*) for t—> 1 if S’ begins with a discontinuity. 
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The problem of minimising the second integral of (12) for given restrictions 
on body geometry is a difficult one. If the method of this note proves worthwhile, 
some progress may be made with this problem by treating the body as a quasi- 
cylinder for purposes of drag minimisation: in this connection it may be significant 
that smooth slender bodies of minimum wave drag are particular cases of quasi- 
cylinders of minimum wave drag. This result is due to Heaslet and Fuller“, whose 
methods may prove helpful in the present context. 


5. Wing-Body-Nacelle Combinations 


To illustrate the application of the foregoing method to configurations with 
engine nacelles, we consider a configuration having a single nacelle whose axis is 
r=ry. The nacelle source strength, gx, is calculated in the same way as that of the 
body, gg. The following notation is now adopted. 


dz, = hwy (z), 


TR) | gx 
(B?r,? - (z-z,)7]J 


dS, + [B'rs? 
|z-2,| <Bry 


(14) 


oo 


Then, under the assumptions of Section 1, the drag is found to be 


D=Dwing + D{ 23+ hwsx} D{hwsx} +D{gex+kw} D{kw). (17) 
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Erratum 
May 1958 (Vol. IX, Part 2) 
In the paper by F. Walkden “ The Shock Pattern of a Wing-Body Combinatio 


Far from the Flight Path” the line immediately after equation (55) on p. 183 i 
misprinted. This line should read 


“In relation (55), s represents the distance along a ray, where the rays sta 
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